如何用R语言对一组样本数据进行拟合求出密度函数?

Python011

如何用R语言对一组样本数据进行拟合求出密度函数?,第1张

曲线拟合:(线性回归方法:lm)

1、x排序

2、求线性回归方程并赋予一个新变量

z=lm(y~x+I(x^2)+...)

3、plot(x,y)#做y对x的散点图

4、lines(x,fitted(z))#添加拟合值对x的散点图并连线

曲线拟合:(nls)

lm是将曲线直线化再做回归,nls是直接拟合曲线。

需要三个条件:曲线方程、数据位置、系数的估计值。

如果曲线方程比较复杂,可以先命名一个自定义函数

例:

f=function(x1, x2, a, b) {a+x1+x2^b}

result=nls(x$y~f(x$x1, x$x2, a, b), data=x, start=list(a=1, b=2))

#x可以是数据框或列表,但不能是矩阵

#对系数的估计要尽量接近真实值,如果相差太远会报错:“奇异梯度”

summary(result) #结果包含对系数的估计和p值

根据估计的系数直接在散点图上使用lines加曲线即可。

曲线拟合:(局部回归)

lowess(x, y=NULL, f = 2/3, iter = 3)

#可以只包含x,也可使用x、y两个变量

#f为窗宽参数,越大越平滑

#iter为迭代次数,越大计算越慢

loess(y~x, data, span=0.75, degree=2)

#data为包含x、y的数据集;span为窗宽参数

#degree默认为二次回归

#该方法计算1000个数据点约占10M内存

举例:

x=seq(0, 10, 0.1)y=sin(x)+rnorm(101)#x的值必须排序

plot(x,y) #做散点图

lines(lowess(x,y)) #利用lowess做回归曲线

lines(x,predict(loess(y~x))) #利用loess做回归曲线,predict是取回归预测值

z=loess(y~x)lines(x, z$fit) #利用loess做回归曲线的另一种做法

R语言实际上是函数的集合,用户可以使用base,stats等包中的基本函数,也可以自己编写函数完成一定的功能。但是初学者往往认为编写R函数十分困难,或者难以理解。这里对如何编写R函数进行简要的介绍。

函数是对一些程序语句的封装。换句话说,编写函数,可以减少人们对重复代码书写,从而让R脚本程序更为简洁,高效。同时也增加了可读性。一个函数往往完成一项特定的功能。例如,求标准差sd,求平均值,求生物多样性指数等。R数据分析,就是依靠调用各种函数来完成的。但是编写函数也不是轻而易举就能完成的,需要首先经过大量的编程训练。特别是对R中数据的类型,逻辑判别、下标、循环等内容有一定了解之后,才好开始编写函数。 对于初学者来说,最好的方法就是研究现有的R函数。因为R程序包都是开源的,所有代码可见。研究现有的R函数能够使编程水平迅速提高。

R函数无需首先声明变量的类型,大部分情况下不需要进行初始化。一个完整的R函数,需要包括函数名称,函数声明,函数参数以及函数体几部分。

函数名称,即要编写的函数名称,这一名称就作为将来调用R函数的依据。

2. 函数声明,包括 <- function, 即声明该对象的类型为函数。

3. 函数参数,这里是输入的数据,函数参数是一个虚拟出来的一个对象。函数参数所等于的数据,就是在函数体内部将要处理的值,或者对应的数据类型。 函数体内部的程序语句进行数据处理,就是对参数的值进行处理 ,这种处理只在调用函数的时候才会发生。函数的参数可以有多种类型。R help的界面对每个函数,及其参数的意义及所需的数据类型都进行了说明。

4. 函数体

常常包括三部分.

(1). 异常处理

输入的数据不能满足函数计算的要求,或者类型不符, 这时候一定要设计相应的机制告诉用户,输入的数据在什么地方有错误。 错误又分为两种。

第一种, 如果输入的数据错误不是很严重,可以经过转换,变为符合处理要求的数据时, 此时只需要给用户一个提醒,告知数据类型不符,但是函数本身已经 进行了相应的转换。

第二种,数据完全不符合要求,这种情况下,就 要终止函数的运行,而告知因为什么,函数不能运行。这样,用户在 使用函数的情况先才不至于茫然。

(2). 运算过程

包括具体的运算步骤。 运算过程和该函数要完成的功能有关。

R运算过程中,应该尽量减少循环的使用,特别是嵌套循环。R提供了 apply,replicate等一系列函数,来代替循环,应该尽量应用这些函数, 提高效率。 如果在R中实在太慢,那么核心部分只能依靠C或者Fortran 等语言编写,然后再用R调用这些编译好的模块,达到更高的效率。

运算过程中,需要大量用到if等条件作为判别的标准。if和while都是需要数据TRUE/FALSE这样的逻辑类型变量,这就意味着,if内部,往往是对条件的判别,例如 is.na, is.matrix, is.numeric等等,或者对大小的比较,如,if(x >0), if(x == 1), if(length(x)== 3)等等。if后面,如果是1行,则花括号可以省略,否则就必须要将所有的语句都放在花括号中。这和循环是一致的。

例子:

## if与条件判断

fun.test <- function(a, b, method = "add"){

if(method == "add") { ## 如果if或者for/while;

res <- a + b ## 等后面的语句只有一行,则无需使用花括号。

}

if(method == "subtract"){

res <- a - b

}

return(res) ## 返回值

}

### 检验结果

fun.test(a = 10, b = 8, method = "add")

fun.test(a = 10, b = 8, method = "substract")

用的最多的,是求均值的mean()函数,当然这里也要提到,像sum()这种求和函数,

还有sd(x) 标准差函数,var(x) 方差函数。min()求最小值,max()求最大值。

我们来具体试试,这里使用一个向量:

test<-c(2,4,5,23,199,25,78,90,12)

求最大值

>max(test)

[1] 19

求最小值

>min(test)

求和

>sum(test)

[1] 43

求标准差,求方差

>sd(test)

[1] 65.01154

>var(test)

[1] 4226.

在来试试最重要的均值

>mean(test)

[1] 48.66667

另外中位数计算。使用median()函数

>median(test)

[1] 23

如果给定一种概率分布,通常会有四类计算问题:

计算其概率密度density (d)计算其概率分布probability(p)计算其百分位数quantile (q)随机数模拟random (r)上面四类计算对应的英文首字母,就是R语言类率分布函数的开头字母。

比如说,正态分布是norm的化,那密度函数就是dnorm(),分布函数就是pnorm(),

更有用的是用相应分布生成随机数,比如rnorm(),就会生成服从正态分布的随机数。

比如我们生成100个服从正态分布的随机数

rnorm(100)

[1] -9.064408e-01 1.026560e+00 -1.097470e+00 1.055395e+00 9.377175e-01

[6] -2.080103e-01 -3.092396e-01 -8.739942e-01 -1.242774e+00 1.102486e+00

[11] 1.082092e+00 -1.695528e+00 -5.930809e-01 -2.100800e-01 8.253859e-01

[16] -1.112551e+00 -3.960474e-01 -9.354820e-01 7.291608e-01 -3.773510e-01

[21] -3.438082e-01 -7.378688e-02 -9.047609e-01 -1.036344e+00 9.485103e-01

[26] -3.437985e-01 -2.145275e-02 1.350098e+00 -1.283633e+00 3.767240e-01

[31] 1.169566e+00 -4.325399e-01 -9.215626e-02 3.839357e-01 3.045491e-01

......

我们再用相应的频率分布直方图来看一下,这些生成的随机数:

hist(rnorm(100))

R就画出了这些随机数的频率分布图