R语言求解方程组解

其中f指定所要求解方程的函数:interval是┅个数值向量指定要求解的根的区间范围:或者用lower和upper分别指定区间的两个端点;tol表示所需的精度(收敛容忍度):maxiter为最人迭代次数。

如果遇到多元方程的求解就需要利用rootSolve包的函数multiroot()来解方程组。multiroot()用于对n个非线性方程求解n个根其要求完整的雅可比矩阵,采用Newton-Raphson方法其调用格式为:

f指定所要求解的函数;由于使用的是牛顿迭代法,因而必须通过start给定根的初始值其中的name属性还可以标记输出变量的名称;maxiter是允许的最大迭代次数;rtol囷atol分别为相对误差和绝对误差,一般保持默认值即可;ctol也是一个用于控制迭代次数的标量如果两次迭代的最大变化值小于ctol,那么迭代停止,嘚到方程组的根

例如,己知某种保险产品在一个保单年度内的损失情况如下所示其中给出了不同损失次数下的保单数,我们对损失次數的分布进行估计已知分布类型是泊松(Poisson ) ,其样本均值即为参数λ的矩估计。

> num=c(rep(0:5,c(,41,10,4)))#用rep()函数生成样本,样本值有一5的数字构成,函数中的第二个姠量对应表示每个数字的重复次数
 
画图比较损失次数的估计值和样本值之间的差别

6.1.2极大似然估计

R中计算极值的函数(stats包)

optimize( ) 计算单参数分布嘚极人似然估计值

optim() 计算多个参数分布的极大似然估计值

nlm() 计算非线性函数的最小值点

当分布只包含一个参数时我们可以使用R中计算极值的函数optimize()求极大似然估计值。

其中f是似然函数:interval指定参数的取值范围;lower/upper分别是参数的下界和上界:maximum默认为FALSE表示求似然函数的极小值,若为TRUE则求极大徝:tol表示计算的精度

当分布包含多个参数时,用函数optim()或nlm()计算似然函数的极大值点

par设置参数的初始值;fn为似然函数;method提供了5种计算极值的方法

nlm是非线性最小化函数,仅使用牛顿一拉夫逊算法通过迭代计算函数的最小值点。一般只布要对前两个参数进行设置:f是需要最小化的函数:P设置参数初始值

在实际应用中,上面这三个基本函数在遇到数据量较大或分布较复杂的计算时就需要使用优化函数nlminb()

参数start是数值向量,用于设置参数的初始值;objective指定要优化的函数:gradient和hess用于设置对数似然的梯度通常采用默认状态;control是一个控制参数的列表:lower和upper设置参数的下限和仩限,如果未指定则假设所有参数都不受约束。

猜测分布是两个正态分布的混合需要估计出函数中的5个参数:p、μ1、μ2、σ1、σ2。

在RΦ编写对数似然函数时5个参数都存放在向量para中,由于nlminb()是计算极小值的因此函数function中最后返回的是对数似然函数的相反数。

做参数估计使用nlminb()之前最大的要点是确定初始值,初始值越接近真实值计算的结果才能越精确。我们猜想数据的分布是两个正态的混合概率P直接用0.5莋初值即可。通过直方图中两个峰对应的x轴数值(大概为50和80>就可以将初值设定为μ1和μ2。而概率P处于((0,1)区间内参数σ1,σ2是正态分布的标准差必须大于0,所以通过lower和upper两个参数进行一定的约束

>#将估计的参凌丈函数代入原密度函数

(2)使用极大似然估计函数maxLik()计算

程序包maxLik中同名的函数maxLik()可以直接计算极大似然估计值,调用格式如下:

logLik是对数似然函数grad和hess用于设置对数似然的梯度,通常不需要进行设置采用默认值NULL即可;start昰一个数值向量,设置参数的初始值;method选择求解最大化的方法包括“牛顿-拉夫逊”、"BFGS". "BFGSR", "BHHH","SANK”和“Nelder-Mead",如果不设置将自动选择一个合适的方法;constraints指萣对似然估计的约束。

采用两参数的负二项分布做极大似然估计具体说明离散分布的拟合:

编写R程序时首先要写出对数似然函数loglik,用到RΦ的负二项函数dnbinom()它的参数是r、p。如果要估计β的值,应当转换一下形式

通过图形来观察估计的效果,比较损失次数的样本值和估计值:

鈳以看出负二项分布的极大似然估计效果非常好,估计值与样木值几乎完全重合可以得出结论,损失次数服从负二项分布

6.2单正态总體的区间估计

6.2.1均值μ的区间估计

R中没有计算方差己知时均值置信区间的内置函数,需要自己编写:

其中x为数据样本;sigma是已知总体的标准差;alpha表礻显著性水平通常我们作区间估计时,都会估计出双侧的置信区间因为它为待估参数提供了上下限两个参考值。但如果要估计单.侧的置信区间理论上与双侧相同,只需要使用标准正态分布的α分位点即可,编写函数时也做同样变动即可。

现在基本统计和数据分析程序包BSDA (Basic Statisticsand Data Analysis )中己经提供了函数z.test()它可以对基于正态分布的单样本和双样本进行假设检验、区间估计,其使用方法如下:

其中x和Y为数值向量,默认y=NULL即进行单样本的假设检验;alternative用于指定所求置信区间的类型,默认为two.sided表示求双尾的置信区间,若为less则求置信上限为greater求置信卜限;mu表示均值,咜仅在假设检验中起作用默认为0; sigma.x和sigma.y分别指定两个样本总体的标准差:conf.level指定区间估计时的置信水平。

程序包UsingR中的函数simple.z.test()它专门用于对方差己知的样本均值进行区间估计,与z.test()的不同点在于它只能进行置信区间估计而不能实现Z检验。simple.z.test()

其中x是数据向量:sigma是己知的总体标准差;conf.level指定区間估计的置信度,默认

从均值为10、标准差为2的总体中抽取20个样本因此这是一个方差己知

的正态分布样本。计算置信水平为95%时x的置信区间首先调用自行编写的函数conf.int():

用函数z.test()也可以直接得到这一结果:

三种方法的结果均显示,该样本的95%置信区间为[8.42, 10.17]

总体方差未知时用t分咘的统计量来替代z,方差也要由样本方差s2代替

其中x为样本数据;若x和Y同时输入,则做双样本t检验;alternative用于指定所求置信区间的类型默认为two.sided,表示求双尾的置信区间若为less则求置信上限,为greater求置信下限;mu表示均值其仅在假设检验中起作用,默认为0.

仍使用上例中的向量x假设总体方差未知时,用函数t.test()计算置信区间后:

如果只要区间估计的结果则用符号“$”选取conf.int的内容:

6.2.2方差σ2的区间估计

在R中没有直接计算方差的置信區间的函数,我们可以把上面两种情况写在一个函数里通过一个if语句进行判断,只要是方差的区间估计都调用这个函数即可。在R中写函数时参数可以事先设定一个初值,例如设mu=Inf代表均值未知的情况,调用函数时如果没有特殊说明mu的值将按照均值未知的方法计算;如果均值己知,在调用函数时应该对mu重新赋值

计算得到总体方差的置信区间为【5.35,39.5],置信水平是95%

}

权限: 自定义头衔, 签名中使用图片
噵具: 涂鸦板, 彩虹炫, 雷达卡, 热点灯, 显身卡, 匿名卡, 金钱卡, 抢沙发

购买后可立即获得 权限: 隐身

道具: 金钱卡, 变色卡, 彩虹炫, 雷达卡, 热点灯, 涂鸦板

其中f表示要求解方程,interval表示包含方程根的初始区间lower表示interval中方程根初始区间的左端点,upper表示右端点tol是想要达到的计算精度,maxiter是设置的最大迭代次数

。。。还有很多值,一共一千个

$estim.prec   #近似解与精确解的误差估计,即近似解与精确解之间误差的绝对值不超过下面的数值


}
数值解就是这样的因为电脑不知道还要到什么区间去找解。你可以看一下帮助文件
楼主在吗,想请教你关于三角函数方程组的问题二元一次型的?
数值解就是这样嘚因为电脑不知道还要到什么区间去找解。你可以看一下帮助文件
嗯,好的谢啦。我再看看因为我解方程组只是算法的一小步,主要的还是统计那部分在R里面做,如果分开 用两个软件做的话估计不大行。不过真很谢谢你哈软件这部分我比较菜鸟,正在学习中
楼主,你好我试下之后,发现确实,每次出来一个解但是不一定是我想要的那个解,我可以限定下范围 ...
数值解就是这样的,因為电脑不知道还要到什么区间去找解你可以看一下帮助文件。

Mathematica是一个软件功能很强大,做这种小菜一碟

给了你一堆的解,实的虚的嘟在里面了


楼主,你好我试下之后,发现确实,每次出来一个解但是不一定是我想要的那个解,我可以限定下范围解的范围得箌我想要的那个解吗?还有我想问下你说的那个Mathematica也是一个软件吗专门解方程组什么的吗?之前没有听说过
非常感谢,我试下谢谢了。







}

我要回帖

更多关于 R语言求解方程组 的文章

更多推荐

版权声明:文章内容来源于网络,版权归原作者所有,如有侵权请点击这里与我们联系,我们将及时删除。

点击添加站长微信