R语言自定义极大似然估计函数、假设检验函数

目录

一、自定义极大似然函数

1、求出似然函数，以正太分布、指数分布为例
正太分布似然函数：

指数分布似然函数（假设x>0）：

2、函数代码：
正太分布求极大似然函数的函数代码

norm.fun<- function(theta,x) # 创建似然函数
{	
	mu<- theta[1]
	sigma<- theta[2] 	# 注意，这里的xigma是方差，不是标准差
	n<- length(x)
	# 似然函数，再求其对数
	logL<- -0.5*n*log(2*pi)-0.5*n*log(sigma)-(1/(2*sigma))*sum((x-mu)**2) 
	return(-logL)  # 因为R中只有最小化函数optim(),我们只需要参数的值，最大化 logL 和最小化 -logL 一致的

theta0=c(0,1) #初始化两个参数作为迭代初始值
X=1:50 #抽出的样本（这里我以1至50作为数据样本）

result=optim(theta0,norm.fun,x=X) 
# 用R中内置的optim()最小化函数对 似然函数优化，求出最优值的两个参数u、sigma
# 参数的值保存在result$'par' 中，有两个值，第一个是u 第二个是sigma（方差）

}

同样指数函数的函数代码：

lamda.fun<- function(theta,x)
{
	n<- length(x)
	logL<- n*log(theta)-theta*sum(x)
	return(-logL)
}

X=1:50 
lamda0<- 1e-2  # 初始化迭代的参数
result1<- optim(lamda0,lamda.fun,x=X)

运行会有一些警告，说optim（）优化的结果不一定正确，但这里是正确的，不影响

二、自定义假设检验函数

基本库中未提供的正态总体方差已知时的均值的区间估计及假设检验函数，就拿这个函数来做例子吧
1、问题分析

**注意，我这里的这里的sigma(那个像6的字母）代表方差，不是标准差（因为之前打公式忘记加一个平方上去，先不管它了）；
单边检验，参数less表示备选假设小于，原假设是大于；参数greater表示备选假设大于，原假设小于；检验区间由枢轴量N推导而出，Inf表示无穷;
当带入数据样本后得到的检验量N属于接受域（拒绝域之外），那么就接受原假设
双边检验的接受域：

单边检验less 接受域：

单边检验 greater 接受域：

…………………………………………………………………………
下面是求置信区间：

** 再次提醒 ，我这里的这里的sigma(那个像6的字母）代表方差，不是标准差，下面的函数代码中的参sigma也是方差

2、函数构造

N.test<- function(x,mu,sigma,a=0.05,alternative='two.side')
#参数依次：数据样本、检验的均值、方差、执行水平a（默认0.05），备选假设是（双边、单边）
{  xbar=mean(x)  # 均值
   n=length(x)	 
   
   if(alternative=='two.side') 
	{
	   N=qnorm(a/2)  # 分位数
	   area1=xbar-(sigma/n)**0.5*N
	   area2=xbar+(sigma/n)**0.5*N
	   area=c(area1,area2)
	   # 上面是由公式写出置信区间
	   
         NT=(xbar-mu)/(sigma/n)**0.5 # 计算N检验值
         
	   TF='reject' # 先默认TF 为拒绝
	  
	   if ( abs(NT)<abs(N))  # 检查 N检验值是否在接受域内
		{
		   TF='accept'  # 如果在接受域内，TF 变为接受 
		}
	result=list('accept_or_reject'=TF,'interval estimation'=area)
	# 返回 是否接受原假设 和 区间估计

	}
# 下面的 less 和greater 思路和two.side 的一样，换公式即可
   if(alternative=='less')
	{
	   N=qnorm(a)
	   area1=xbar+(sigma/n)**0.5*N
	   area2=Inf  # 无穷
	   area=c(area1,area2)
	   
         NT=(xbar-mu)/(sigma/n)**0.5
	   TF='reject'
	   if (NT>N)
		{
		   TF='accept'
		}
	result=list('accept_or_reject'=TF,'interval estimation'=area)

	}

   if(alternative=='greater')
	{
	   N=-qnorm(a)   #N应该是a的上分位，对下分位求负数是相等的
	   area1=-Inf
	   area2=xbar-(sigma/n)**0.5*N

	   area=c(area1,area2)
	   
         NT=(xbar-mu)/(sigma/n)**0.5
	   TF='reject'
	   if (NT<N)
		{
		   TF='accept'
		}
	result=list('accept_or_reject'=TF,'interval estimation'=area)
	}
	
   return(result)
}

3、运行结果，检验函数：
检验H0: mu=0,用双边

结果是接受的，符合实际

检验H0: mu>1，此时备选假设是mu<1,用单边检验，参数less

结果是拒绝H0，即不认为mu>1，这也是符合实际的

检验H0： mu<1,此时备选假设是mu>1，参数greater

结果是接受H0的，认为mu<1，这也是符合实际的

后面陆续做了多组数据的检验，检验效果都合理。