R 的 integrate 問題

December 21, 2010

測試一下:
\[ \Phi(\color{red}x) = \int_{-\infty}^{\color{red}x} \phi(z)dz\]

> integrate(dnorm,0,4384)
0.5 with absolute error < 5.2e-05
> integrate(dnorm,0,4385)
6.124136e-05 with absolute error < 0.00012

理論上應該是要遞增才對﹐ 一下掉到快變零,這是為什麼?

Advertisements

MathML stretch symbol Test

October 25, 2010

E
(

X
k

)
=

R

1

Γ

ν
+
1

2

Γ

ν
+
1

2

π
ν

x
k

ν

ν
+

x
2

ν
+
1

2

d
x

Nothing selected

double gsl_cdf_noncentral_F_P(double x, double nu1, double nu2, double ncp) {
  double resultP = gsl_sf_exp(-ncp/2.0);
  double mediSum = 0.0;
  double mediSum1 = 1.0;
  int i = 0;
  while (fabs(mediSum1 - mediSum) > 0) {
    mediSum1 = mediSum;
    mediSum = mediSum + pow(ncp/2.0, (double) i)*gsl_cdf_beta_Q (nu2/(nu2+nu1*x), nu2/2.0, nu1/2.0+i )/gsl_sf_fact(i);
    i = i + 1;
  }
  resultP = resultP*mediSum;
  return resultP;
}

不廢話,code:

/*
** testing.c
** 
** Made by (Gong-Yi Liao)
** Login   
** 
** Started on  Sun Mar 14 22:12:22 2010 Gong-Yi Liao
** Last update Sun May 12 01:17:25 2002 Speed Blue
*/

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

double gsl_cdf_noncentral_chisq_P(double x, double nu, double ncp) {
  double resultP = gsl_sf_exp(-ncp/2.0);
  double mediSum = 0.0;
  double mediSum1 = 1.0;
  int i = 0;
  while (fabs(mediSum1 - mediSum) > 0) {
    mediSum1 = mediSum;
    mediSum = mediSum + pow(ncp/2.0, (double) i)*gsl_sf_gamma_inc_P (i+nu/2.0, x/2.0)/gsl_sf_fact(i);
    i = i + 1;
  }
  resultP = resultP*mediSum;
  return resultP;
}

int main(int argc, char *argv[]) {
  printf("%e\n", gsl_cdf_fdist_P (5.2, 7.4, 5.8));
  printf("%e\n", gsl_cdf_chisq_P (5.2, 7.4));
  printf("%e\n", gsl_cdf_noncentral_chisq_P (5.2, 7.4, 14.1));
  return 0;
}

讀書讀累了, 把之前的gsl_clapack_dgesdd_() 作一點小修改, 把$\mathbf U$, $\mathbf s$, 和$\mathbf V’$ 的變數設定的責任丢回給使用者
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <f2c.h>
#include <clapack.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

int gsl_clapack_dgesdd_(gsl_matrix *mtrxA, gsl_matrix *mtrxU, gsl_vector *vecS, gsl_matrix *mtrxV) {
  int works = 0;
  int M = mtrxA->size1;
  int N = mtrxA->size2;
  double *A = (double*) malloc(sizeof(double)*M*N);
  memcpy(A, mtrxA->data, N*M);
  char jobz = 'S';
  int lda = max(1, M);
  int ldu = max(M, N);
  int ldvt = min(M, N);
  double *s = vecS->data;
  int lwork = 4*min(M,N) + max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N));
  int *iwork = malloc(sizeof(int)*8*max(M, N));
  double *wk = malloc(sizeof(double)*lwork);
  double *uu = mtrxU->data;
  double *vt = mtrxV->data;
  int info;

  if (dgesdd_(&jobz, &M, &N, A, &lda, s, uu, &ldu, vt, &ldvt, wk, &lwork, iwork, &info) == 1) 
    works = 1;

  free(wk);
  free(iwork);
  free(A);
  
  return works;
}


int main(int argc, char **argv) {
  int M = atoi(argv[1]);
  int N = atoi(argv[2]); 
  gsl_rng *rng;
  const gsl_rng_type *rngType;    
  gsl_rng_env_setup();
  rngType = gsl_rng_default;
  rng = gsl_rng_alloc(rngType);
  
  gsl_matrix *A = gsl_matrix_alloc(M, N);

  int i = 0;
  int j = 0;
  for (i = 0; i < M; i++) 
    for (j = 0; j < N; j++) 
      gsl_matrix_set(A, i, j, gsl_ran_ugaussian(rng)); 

  gsl_matrix *U = gsl_matrix_alloc(max(M,N), min(M,N));
  gsl_vector *s = gsl_vector_alloc(min(M, N));
  gsl_matrix *V = gsl_matrix_alloc(min(M, N), N);

  if (!(gsl_clapack_dgesdd_(A, U, s, V) == 0))
    printf("Error!\n"); 

  printf("%e\n", gsl_matrix_get(U, 20, 20));

  gsl_rng_free(rng);

  return 0;
}

dgesdd_ from CLAPACK:


#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <f2c.h>
#include <clapack.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

typedef struct  {
  gsl_matrix *svdUmatrix;
  gsl_matrix *svdVmatrix;
  gsl_vector *svdSvector;
} gsl_lapack_svd_ds;

gsl_lapack_svd_ds* gsl_clapack_dgesdd_(gsl_matrix *mtrxA) {
  int M = mtrxA->size1;
  int N = mtrxA->size2;
  double *A = malloc(sizeof(double)*M*N);
  memcpy(A, mtrxA->data, N*M);
  char jobz = 'S';
  int lda = max(1, M);
  int ldu = max(M, N);
  int ldvt = min(M, N);
  double *s = malloc(sizeof(double)*min(M, N));
  int lwork = 4*min(M,N) + max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N));
  int *iwork = malloc(sizeof(int)*8*max(M, N));
  double *wk = malloc(sizeof(double)*lwork);
  double *uu = malloc(sizeof(double)*ldu*min(M, N));
  double *vt = malloc(sizeof(double)*ldvt*N);
  int info;

  dgesdd_(&jobz, &M, &N, A, &lda, s, uu, &ldu, vt, &ldvt, wk, &lwork, iwork, &info);

  gsl_lapack_svd_ds *gsl_svd_sol1 = (gsl_lapack_svd_ds*) malloc(sizeof(gsl_lapack_svd_ds));
  
  gsl_matrix_view svdUview = gsl_matrix_view_array(uu, ldu, min(M, N));
  gsl_svd_sol1->svdUmatrix = gsl_matrix_alloc(ldu, min(M,N));
  gsl_matrix_memcpy(gsl_svd_sol1->svdUmatrix, &svdUview.matrix);

  gsl_matrix_view svdVview = gsl_matrix_view_array(vt, ldvt, N);
  gsl_svd_sol1->svdVmatrix = gsl_matrix_alloc(ldvt ,N);
  gsl_matrix_memcpy(gsl_svd_sol1->svdVmatrix, &svdVview.matrix);

  gsl_vector_view svdSview = gsl_vector_view_array(s, min(M, N));
  gsl_svd_sol1->svdSvector = gsl_vector_alloc(min(M, N));
  gsl_vector_memcpy(gsl_svd_sol1->svdSvector, &svdSview.vector);

  free(wk);
  free(iwork);
  free(A);
  
  return gsl_svd_sol1;
}


int main(int argc, char **argv) {
  int M = atoi(argv[1]);
  int N = atoi(argv[2]); 
  gsl_rng *rng;
  const gsl_rng_type *rngType;    
  gsl_rng_env_setup();
  rngType = gsl_rng_default;
  rng = gsl_rng_alloc(rngType);
  
  gsl_matrix *A = gsl_matrix_alloc(M, N);

  int i = 0;
  int j = 0;
  for (i = 0; i < M; i++) 
    for (j = 0; j svdVmatrix, 20, 20));

  gsl_rng_free(rng);

  return 0;
}