柯西分布正割法 C++ Rcpp MLE/Root

Question

我对 C++ 和 Rcpp 集成还很陌生。我需要创建一个使用 C++ 和 R 集成的程序来找到 Cauchy 分布的 MLE/root。

以下是我到目前为止的代码。

#include <Rcpp.h> 
#include <math.h>
#include <iostream>
#include <cstdlib>
using namespace std;
using namespace Rcpp;

// [[Rcpp::export]]
double Cauchy(double x, double y); //Declare Function
double Cauchy(double x,double y)   //Define Function
{
    return 1/(M_PI*(1+(pow(x-y,2)))); //write the equation whose roots are to be determined x=chosen y=theta 
}

using namespace std;
// [[Rcpp::export]]

int Secant (NumericVector x){
   NumericVector xvector(x) ; //Input of x vector
   double eplison= 0.001; //Threshold
   double a= xvector[3]; //Select starting point
   double b= xvector[4];//Select end point
   double c= 0.0; //initial value for c
   double Theta= 10.6; //median value for theta estimate
   int noofIter= 0; //Iterations
   double error = 0.0; 

   if (std::abs(Cauchy(a, Theta) <(std::abs(Cauchy(a, Theta))))) {  
      do {
         a=b;
         b=c; 
         error= (b-(Cauchy(b, Theta)))*((a-b)/(Cauchy(a, Theta)-Cauchy(b, Theta))); 
         error= Cauchy(c,Theta); 
         //return number of iterations
         noofIter++;

      for (int i = 0; i < noofIter; i += 1) {
         cout << "The Value is " << c << endl;
         cout << "The Value is " << a << endl;
         cout << "The Value is " << b << endl;
         cout << "The Value is " << Theta << endl;
       }
      } while (std::abs(error)>eplison);
    }

    cout<<"\nThe root of the equation is occurs at "<<c<<endl;    //print the root
    cout << "The number of iterations is " << noofIter;
    return 0;
}

通过一些修改，程序要么进入永无止境的循环，要么 returns 一个无限小的值。

我对这个数学的理解是有限的。因此，我们将不胜感激任何帮助或更正。

我们作为输出给出的 X 向量是

x <- c( 11.262307 , 10.281078 , 10.287090 , 12.734039 ,
         11.731881 , 8.861998 , 12.246509 , 11.244818 ,
         9.696278 , 11.557572 , 11.112531 , 10.550190 ,
         9.018438 , 10.704774 , 9.515617 , 10.003247 ,
         10.278352 , 9.709630 , 10.963905 , 17.314814)

使用之前的 R 代码我们知道这个分布的 MLE/Root 出现在 10.5935 大约

用于获取此 MLE 的代码是

 optimize(function(theta)-sum(dcauchy(x, location=theta, 
     log=TRUE)),  c(-100,100))

谢谢！

Answer 1

使用 optimize() 函数，您可以直接搜索可能性的极值。另一种方法是将求根算法（例如割线法）与（对数）似然的导数一起使用。从 Wikipedia 我们得到了我们必须解决的公式。在 R 中，这可能是这样的：

x <- c( 11.262307 , 10.281078 , 10.287090 , 12.734039 ,
        11.731881 , 8.861998 , 12.246509 , 11.244818 ,
        9.696278 , 11.557572 , 11.112531 , 10.550190 ,
        9.018438 , 10.704774 , 9.515617 , 10.003247 ,
        10.278352 , 9.709630 , 10.963905 , 17.314814)

ld <- function(sample, theta){
  xp <- outer(sample, theta, FUN = "-")
  colSums(xp/(1+xp^2))
}
uniroot(ld, sample = x, lower = 0, upper = 20)$root
#> [1] 10.59724

请注意，对数似然的导数在两个参数上都被向量化了。这样可以轻松绘图：

theta <- seq(0, 20, length=500)
plot(theta, ld(x, theta), type="l",
     xlab=expression(theta), ylab=expression(ld(x, theta)))

从该图中我们已经看到，要找到正确的正割法工作起点将很棘手。

让我们将其移至 C++（准确地说是 C++11）：

#include <Rcpp.h>
// [[Rcpp::plugins(cpp11)]]

Rcpp::List secant(const std::function<double(double)>& f, 
              double a, double b, int maxIterations, double epsilon) {
  double c(0.0);
  do {
    c = b * (1 - (1 - a/b) / (1 - f(a)/f(b)));
    a = b;
    b = c;
  } while (maxIterations-- > 0 && std::abs(a - b) > epsilon);
  return Rcpp::List::create(Rcpp::Named("root") = c,
                            Rcpp::Named("f.root") = f(c),
                            Rcpp::Named("converged") = (maxIterations > 0));
}

// [[Rcpp::export]]
Rcpp::List mleCauchy(const Rcpp::NumericVector& sample, double a, double b,
                     int maxIterations = 100, double epsilon = 0.0001) {
  auto f = [&sample](double theta) {
    Rcpp::NumericVector xp = sample - theta;
    xp = xp / (1 + xp * xp);
    return Rcpp::sum(xp);
  };
  return secant(f, a, b, maxIterations, epsilon);
}


/*** R
x <- c( 11.262307 , 10.281078 , 10.287090 , 12.734039 ,
        11.731881 , 8.861998 , 12.246509 , 11.244818 ,
        9.696278 , 11.557572 , 11.112531 , 10.550190 ,
        9.018438 , 10.704774 , 9.515617 , 10.003247 ,
        10.278352 , 9.709630 , 10.963905 , 17.314814)

mleCauchy(x, 11, 15)
#-> does not converge
mleCauchy(x, 11, 14)
#-> 10.59721
mleCauchy(x, mean(x), median(x))
#-> 10.59721
*/

secant() 函数适用于任何 std::function that takes a double as argument and returns a double. Such a function is than defined as lambda function 取决于提供的示例值。正如预期的那样，只有当从接近正确值的值开始时，才能获得正确的根。

Lambda 函数乍一看可能有点混乱，但它们非常接近我们在 R 中使用的函数。这里用 R 编写相同的算法：

secant <- function(f, a, b, maxIterations, epsilon) {
  for (i in seq.int(maxIterations)) {
    c <- b * (1 - (1 - a/b) / (1 - f(a)/f(b)))
    a <- b
    b <- c
    if (abs(a - b) <= epsilon)
      break
  }
  list(root = c, f.root = f(c), converged = (i < maxIterations))
}

mleCauchy <- function(sample, a, b, maxIterations = 100L, epsilon = 0.001) {
  f <- function(theta) {
    xp <- sample - theta
    sum(xp/(1 + xp^2))
  }
  secant(f, a, b, maxIterations, epsilon)
}

x <- c( 11.262307 , 10.281078 , 10.287090 , 12.734039 ,
        11.731881 , 8.861998 , 12.246509 , 11.244818 ,
        9.696278 , 11.557572 , 11.112531 , 10.550190 ,
        9.018438 , 10.704774 , 9.515617 , 10.003247 ,
        10.278352 , 9.709630 , 10.963905 , 17.314814)
mleCauchy(x, 11, 12)
#-> 10.59721

R 函数 f 和 lambda 函数 f 从定义它们的环境中获取向量 sample。在 R 中，这是隐式发生的，而在 C++ 中，我们必须明确告知应该捕获该值。数字 theta 是调用函数时提供的参数，即以 a 和 b.

开头的根的连续估计

柯西分布正割法 C++ Rcpp MLE/Root

Cauchy Distribution Secant Method C++ Rcpp MLE/Root

c++

math

r

rcpp