R语言求解非线性不等式组

问题见图片附件，急求 解出该非线性不等式组 的R语言程序！！
【初步求解思路：首先选定个较小的n (如n=2)，利用二分法查找出第一个公式在等式成立时C0的值:再将上述值代入第二个式子，以验证所求得的C0是否满足第二个式子; 如不满足，则样本量n的值加 1，直到第一个式子解出的n和C0也满足第二个式子】
【初步程序（报错）：
f_C1<-function(t){
  (pchisq((n-1)*((3*C1+1)*sqrt(n)-t)^2/(9*n*C0^2),df=n-1))*(dnorm(t+sqrt(n))+dnorm(t-sqrt(n)))}
f_C2<-function(t){
  (pchisq((n-1)*((3*C2+1)*sqrt(n)-t)^2/(9*n*C0^2),df=n-1))*(dnorm(t+sqrt(n))+dnorm(t-sqrt(n)))}
C1 <- 4/3; C2 <- 1; alpha <- 0.05; n <- 2
#二分法求解方程的根
S <- function(C0){integrate(f_C1,lower=0, upper=(3*C1+1)*sqrt(n))$value-(1-alpha)}
x0 <- 0.5; x1 <- 2
while(abs(x0-x1)>1e-5){
  if(S(x0)*S(x1)>0){
    print("no root")
    break
  }
  x2 = (x0+x1)/2
  if(S(x2)==0){cat(x2);break}
  if(S(x0)*S(x2)<0){
    x1 = x2
  }else{
    x0 = x2
  }
}
##此时报错：Error in pchisq((n - 1) * ((3 * C1 + 1) * sqrt(n) - t)^2/(9 * n *  : object 'C0' not found                                                                                                                  】

在程序中，C0是待求解的未知量，没有事先给定值，因此无法在函数中直接使用。可以通过将C0作为参数传递给函数来解决这个问题。另外，为了避免后面的函数定义中对n和C0的使用报错，需要将这两个变量的定义放到前面。

1. C1 <- 4/3
   
2. C2 <- 1
   
3. alpha <- 0.05
   
4. n <- 2
   
5. 
   
6. f_C1 <- function(t, n, C0, C1) {
   
7.   (pchisq((n - 1) * ((3 * C1 + 1) * sqrt(n) - t)^2 / (9 * n * C0^2), df = n - 1)) *
   
8.     (dnorm(t + sqrt(n)) + dnorm(t - sqrt(n)))
   
9. }
   
10. 
    
11. f_C2 <- function(t, n, C0, C2) {
    
12.   (pchisq((n - 1) * ((3 * C2 + 1) * sqrt(n) - t)^2 / (9 * n * C0^2), df = n - 1)) *
    
13.     (dnorm(t + sqrt(n)) + dnorm(t - sqrt(n)))
    
14. }
    
15. 
    
16. # 二分法求解方程的根
    
17. S <- function(C0, n, C1, alpha) {
    
18.   integrate(f_C1, lower = 0, upper = (3 * C1 + 1) * sqrt(n), n = n, C0 = C0, C1 = C1)$value - (1 - alpha)
    
19. }
    
20. 
    
21. x0 <- 0.5
    
22. x1 <- 2
    
23. 
    
24. while (abs(x0 - x1) > 1e-5) {
    
25.   if (S(x0, n, C1, alpha) * S(x1, n, C1, alpha) > 0) {
    
26.     print("no root")
    
27.     break
    
28.   }
    
29.   x2 <- (x0 + x1) / 2
    
30.   if (S(x2, n, C1, alpha) == 0) {
    
31.     cat(x2)
    
32.     break
    
33.   }
    
34.   if (S(x0, n, C1, alpha) * S(x2, n, C1, alpha) < 0) {
    
35.     x1 <- x2
    
36.   } else {
    
37.     x0 <- x2
    
38.   }
    
39. }

复制代码

oliyiyi 发表于 2023-2-24 08:59

感谢大神解惑！
能否请大神再帮忙写出循环程序，如【初步求解思路】所述，解出能同时满足两个公式的n和C0？
还有就是，不熟悉页面，不小心点踩了，抱歉。。。。