假设我有从1951年到1980年的100个站点的年降水量数据.在一些论文中,我发现人们将PCA应用于时间序列,然后绘制空间载荷图(值为-1到1),并绘制时间序列个人电脑 例如,https://publicaciones.unirioja.es/ojs/index.php/cig/article/view/2931/2696中的图6 是PC的空间分布.

我prcomp在R中使用函数,我想知道我怎么能做同样的事情.换句话说,如何从prcomp函数的结果中提取"空间模式"和"时间模式" ？谢谢.

set.seed(1234)
rainfall = sample(x=100:1000,size = 100*30,replace = T)
rainfall=matrix(rainfall,nrow=100)
colnames(rainfall)=1951:1980
PCA = prcomp(rainfall,retx=T)

或者,有真实的数据在https://1drv.ms/u/s!AnVl_zW00EHegxAprS4s7PDaYQVr

我正在尝试使用 R 实现多目标遗传算法 (EPR-MOGA) 的进化多项式回归。鉴于此主题没有可用的 R 包，我想知道起点是什么？谢谢。

以下是一些讨论该方法的参考文献。

https://iwaponline.com/jh/article/11/3-4/225/40/Advances-in-data-driven-analysiss-and-modelling

https://iwaponline.com/jh/article/8/3/207/31275/A-symbolic-data-driven-technique-based-on

https://ojs.library.queensu.ca/index.php/wdsa-ccw/article/view/12093

我现在正在使用该extRemes软件包构建一个非固定的GP模型,我发现很难提取参数.

非平稳尺度参数

library(extRemes)
data(Fort)
fit1 <- fevd(Prec, Fort, threshold=0.395,
        scale.fun=~sin(2 * pi * (year - 1900)/365.25) + cos(2 * pi * (year - 1900)/365.25),
        type="GP", use.phi=TRUE, verbose=TRUE)

根据fevd帮助页面,log(scale(y)) = phi(y) = phi0 + phi1 * g1(y) + phi2 * g2(y) + ...

现在,我们有phi0,phi1,并phi2从结果来看,但什么是g1(y)和g2(y)上面的功能？

此外,我们怎么能理解scale.fun的fit1？什么scale.fun=~sin(2 * pi * (year - 1900)/365.25) + cos(2 * pi * (year - 1900)/365.25)代表什么？例如,如果我们使用scale.fun=~Fort$year …

我现在正在R中研究copula，不知道如何找到R中的联合累积分布？

D = c(1,3,2,2,8,2,1,3,1,1,3,3,1,1,2,1,2,1,1,3,4,1,1,3,1,1,2,1,3,7,1,4,6,1,2,1,1,3,1,2,2,3,4,1,1,1,1,2,2,12,1,1,2,1,1,1,3,4)
S = c(1.42,5.15,2.52,2.29,12.36,2.82,1.49,3.53,1.17,1.03,4.03,5.26,1.65,1.41,3.75,1.09,3.44,1.36,1.19,4.76,5.58,1.23,2.29,7.71,1.12,1.26,2.78,1.13,3.87,15.43,1.19,4.95,7.69,1.17,3.27,1.44,1.05,3.94,1.58,2.29,2.73,3.75,6.80,1.16,1.01,1.00,1.02,2.32,2.86,22.90,1.42,1.10,2.78,1.23,1.61,1.33,3.53,10.44)

经过一番探索，我发现Gamma分布是描述上述数据的最佳方法。

library(fitdistrplus)
fg_d <- fitdist(data = Dur, distr = "gamma", method = "mle")
fg_s <- fitdist(data = Sev, distr = "gamma", method = "mle")

然后，我尝试使用VineCopulapackge 选择copula系列：

mydata <- cbind(D=D, S=S)
u1 <- pobs(mydata[,1]) 
u2 <- pobs(mydata[,2])
fitCopula <- BiCopSelect(u1, u2, familyset=NA)
summary(fitCopula) 

结果表明“生存克莱顿”。然后，我尝试构建以下copula：

library(copula)
cop_model <- surClaytonCopula(param = 5.79)

现在，根据下面的方程式（假定E（L）为常数）：

我需要找到给定D和S值的FD（d），FS（s）和C（FD（d），FS（s））。

例如，如果我们采用D = 3和S = 2，则必须找到F（D <= 3），F（S <= 2）和C（D <= 3 andS <= 2）。我想知道如何在R中使用软件包copula吗？

另外，我们如何找到C（D <= 3 …

我现在正在处理一个大型数据集，某些功能可能需要几个小时才能处理。我想知道如何通过进度条或数字（1,2,3，...，100）显示代码的进度。这是一个例子。谢谢。

require(Kendall)
mydata=matrix(rnorm(6000*300),ncol = 300)
result=as.data.frame(matrix(nrow = 6000,ncol = 2))
for (i in 1:6000) {
  abc=MannKendall(mydata[i,])
  result[i,1]=abc$tau
  result[i,2]=abc$sl
}

顺便说一句，我发现链接https://ryouready.wordpress.com/2009/03/16/r-monitor-function-progress-with-a-progress-bar/非常有用。但是，我不知道如何将链接中的代码与自己的函数结合在一起。有人有主意吗？谢谢。这是上面链接的代码。

total <- 20
# create progress bar
pb <- txtProgressBar(min = 0, max = total, style = 3)
for(i in 1:total){
   Sys.sleep(0.1)
   # update progress bar
   setTxtProgressBar(pb, i)
}
close(pb)

我现在正在处理大型数据集,某些功能可能需要数小时才能处理.我想知道如何通过进度条或数字(1,2,3,...,100)显示代码的进度.我想将结果存储为具有两列的数据框.这是一个例子.谢谢.

require(foreach)
require(doParallel)
require(Kendall)

cores=detectCores()
cl <- makeCluster(cores-1)
registerDoParallel(cl)

mydata=matrix(rnorm(8000*500),ncol = 500)
result=as.data.frame(matrix(nrow = 8000,ncol = 2))
pb <- txtProgressBar(min = 1, max = 8000, style = 3)

foreach(i=1:8000,.packages = "Kendall",.combine = rbind) %dopar%         
{
  abc=MannKendall(mydata[i,])
  result[i,1]=abc$tau
  result[i,2]=abc$sl
  setTxtProgressBar(pb, i)
}
close(pb)
stopCluster(cl)

但是,当我运行代码时,我没有看到任何进度条显示,结果不正确.有什么建议吗？谢谢.

我想使用ggplot2在投影地图上绘制一些插值数据,我已经在这个问题上工作了几个星期.希望有人可以帮助我,非常感谢.shapefile和数据可以在https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip?dl=0和https://www.dropbox.com/s/9czvb35vsyf3t28/Mydata.rdata?dl=找到.0.

首先,shapefile最初使用"lon-lat"投影,我需要将其转换为Albers Equal Area(aea)投影.

library(automap)
library(ggplot2)
library(rgdal)
load("Mydata.rdata",.GlobalEnv)
canada2<-readOGR("gpr_000b11a_e.shp", layer="gpr_000b11a_e")
g <- spTransform(canada2, CRS("+proj=aea +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
Borders=ggplot() +geom_polygon(data=g,aes(x=long,y=lat,group=group),fill='white',color = "black")
Borders

我们可以看到,我们可以正确地绘制国家.然后我想用Kriging方法插入数据,代码取自Smoothing out ggplot2 map.

coordinates(Mydata)<-~longitude+latitude
proj4string(Mydata)<-CRS("+proj=longlat +datum=NAD83")
sp_mydata<-spTransform(Mydata,CRS(proj4string(g)))
Krig=autoKrige(APPT~1,sp_mydata)
interp_data = as.data.frame(Krig$krige_output)
colnames(interp_data) = c("latitude","longitude","APPT_pred","APPT_var","APPT_stdev")
interp_data=interp_data[,1:3]
ggplot(data=interp_data, aes(x=longitude, y=latitude)) + geom_tile(aes(fill=APPT_pred),color=NA)

然后我们可以看到插值数据图.

最后我想结合这两个数字然后我得到以下错误: Error: Don't know how to add o to a plot

ggplot(data=interp_data, aes(x=longitude, y=latitude)) + geom_tile(aes(fill=APPT_pred),color=NA)+Borders 

我的问题是:如何在地图上绘制插值数据,然后添加网格线(经度和纬度).另外,我想知道如何剪切插值数据图以适应整个加拿大地图.谢谢您的帮助.

我有一个 netcdf 文件（全局域），我不知道它的投影信息，我想根据 shapefile 提取数据（带有 lon 和 lat）。

目标是找到由 shapefile 表示的域的数据（原始 netcdf 包含全局数据）。此外，最终的数据格式应该是一个数据帧，其中包含 lon、lat 和 data。我认为extractandmask函数会很有用。

netcdf 和 shapefile 可以在https://www.dropbox.com/s/ju96eitzzd0xxg8/precip.2000.nc?dl=0和https://www.dropbox.com/s/8wfgf8207dbh79r/gpr_000b11a_e.zip? dl=0。谢谢你的帮助。

library(rgdal)
library(ncdf4)
library(raster)

shp = readOGR("gpr_000b11a_e.shp")
pre1 = nc_open("precip.2000.nc")
pre1

File precip.2000.nc (NC_FORMAT_NETCDF4_CLASSIC):
 1 variables (excluding dimension variables):
    float precip[lon,lat,time]   
        missing_value: -9.96920996838687e+36
        var_desc: Precipitation
        level_desc: Surface
        statistic: Total
        parent_stat: Other
        long_name: Daily total of precipitation
        cell_methods: time: sum
        avg_period: 0000-00-01 00:00:00
        actual_range: 0
         actual_range: 482.555358886719
        units: mm
        valid_range: 0
         valid_range: 1000 …

我有两个向量，我试图找到向量 1 中的 3 个元素和向量 2 中的 2 个元素的所有唯一组合。我已经尝试过以下代码。

V1 = combn(1:5, 3)   # 10 combinations in total
V2 = combn(6:11, 2)  # 15 combinations in total

如何将V1和V2组合起来，使得总共有10 * 15 = 150种组合？谢谢。

我现在正在使用绘制加拿大地图ggplot2。因为默认的投影方法是“ aea”（阿尔伯斯等面积），所以经度和纬度是地图上的直线。我想知道如何在地图上以“ 110W，100W，90W”和“ 50N，60N，70N”的形式显示经度和纬度。它们应该是曲线。非常感谢。

可从https://www.arcgis.com/home/item.html?id=dcbcdf86939548af81efbd2d732336db下载arcgis shapfile

library(ggplot2)
library(rgdal)
countries<-readOGR("Canada.shp", layer="Canada")
ggplot()+geom_polygon(data=countries,aes(x=long,y=lat,group=group),fill='white',color = "black")

最终结果应该是这样的。