Integral, Measure and Martingale

2020年6月6日星期六

SAS Macro functions – %eval() and %sysevalf()

The %EVAL function evaluates integer arithmetic or logical expressions. %EVAL operates by converting its argument from a character value to a numeric or logical expression. Then, it performs the evaluation. Finally, %EVAL converts the result back to a character value and returns that value.

Operands that contain a period character cause an error when they are part of an integer arithmetic expression. The following examples show correct and incorrect usage, respectively:

%let d=%eval(10+20);     /* Correct usage */

%let d=%eval(10.0+20.0); /* Incorrect usage */

When %EVAL encounters a value containing a period, it displays an error message about finding a character operand where a numeric operand is required.

%let a=1+2;

%let b=10*3;

%let c=5/3;

%let eval_a=%eval(&a);

%let eval_b=%eval(&b);

%let eval_c=%eval(&c);

%put &a is &eval_a;

%put &b is &eval_b;

%put &c is &eval_c;

The %SYSEVALF function evaluates arithmetic and logical expressions using floating-point arithmetic and returns a value that is formatted using the BEST32. format. The result of the evaluation is always text. % SYSEVALF is the only macro function that can evaluate logical expressions that contain floating point or missing values. Specifying a conversion type can prevent problems when %SYSEVALF returns missing or floating point values to macro expressions or macro variables that are used in other macro expressions that require an integer value.

Syntax: %SYSEVALF(expression<, conversion-type>)

Taking an example…
view source
print?

%let a=100;

%let b=101.590;

%let x=%sysevalf(&a+&b);

%let y=%sysevalf(&a +&b, boolean);

%let z=%sysevalf(&a +&b, ceil);

%let v=%sysevalf(&a +&b, floor);

%let w=%sysevalf(&a +&b, int);

%put The result with SYSEVALF is: &x;

%put  The BOOLEAN value is: &y;

%put  The CEIL value is: &z;

%put  The FLOOR value is: &v;

%put  The INTEGER value is: &w;

http://technico.qnownow.com/sas-macro-functions-eval-and-sysevalf-2/

2020年6月2日星期二

parmbuff (or Pbuff) in SAS macro

%macro demo(a=1, b=2)/parmbuff;
%put buffer holds |&syspbuff|;
%put &=a;
%put &=b;
%mend demo;
%demo(a=aa)

%demo(a=silly, d=unknown)

2020年3月13日星期五

Inverse Probability Censoring Weighting Estimation

https://documentation.sas.com/?docsetId=statug&docsetTarget=statug_rmstreg_details13.htm&docsetVersion=15.1&locale=en

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5379776/

2019年11月9日星期六

BLUP for mixed model

http://rpubs.com/enwuliu/548305

This is R code for the example from the paper "That Blup Is a Good Thing: The Estimation of Random Effects.".

Robinson, George K. "That Blup Is a Good Thing: The Estimation of Random Effects." Statistical Science 6, no. 1 (1991): 15-32.

In the book Generalized Additive Models An Introduction with R, the second version, page 81 the function llm that solved the linear mixed model also used Henderson's equations but did not mention the equations.

mixed_model_blup.R
enliu
2019-11-10
Y<-matrix(c(110,100,110,100,100,110,110,100,100),nrow=9,ncol=1)
X<-matrix(c(1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1,1),nrow=9,ncol=3)

Z<-matrix(c(1,0,0,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,0,
0,0,0,0,0,1,1,0,0,
0,1,0,1,1,0,0,1,1),nrow=9,ncol=4)
R<-diag(1,nrow=9,ncol=9)
G<-diag(0.1,nrow=4,ncol=4)
A<-t(X)%*%solve(R)%*%X
B<-t(X)%*%solve(R)%*%Z
C<-t(Z)%*%solve(R)%*%X
D<-t(Z)%*%solve(R)%*%Z+solve(G)

E<-cbind(A,B)
F<-cbind(C,D)

M<-rbind(E,F)

t(X)%*%solve(R)%*%Y
## [,1]
## [1,] 210
## [2,] 310
## [3,] 420
t(Z)%*%solve(R)%*%Y
## [,1]
## [1,] 110
## [2,] 110
## [3,] 220
## [4,] 500
N<-rbind(t(X)%*%solve(R)%*%Y,t(Z)%*%solve(R)%*%Y)

solve(M,N)
## [,1]
## [1,] 105.6386574
## [2,] 104.2757378
## [3,] 105.4584366
## [4,] 0.3964857
## [5,] 0.5203875
## [6,] 0.7569272
## [7,] -1.6738004

/*********Another method from the Generalized Additive Models An Introduction with R, the second version, p81

llm <- function(theta,X,Z,y) {
## untransform parameters...
sigma.b <- exp(theta[1])
sigma <- exp(theta[2])
## extract dimensions...
n <- length(y); pr <- ncol(Z); pf <- ncol(X)
## obtain \hat \beta, \hat b...
X1 <- cbind(X,Z)
ipsi <- c(rep(0,pf),rep(1/sigma.b^2,pr))

b1 <- solve(crossprod(X1)/sigma^2+diag(ipsi),t(X1)%*%y/sigma^2)

## compute log|Z'Z/sigma^2 + I/sigma.b^2|...

ldet <- sum(log(diag(chol(crossprod(Z)/sigma^2 +
diag(ipsi[-(1:pf)])))))
## compute log profile likelihood...
l <- (-sum((y-X1%*%b1)^2)/sigma^2 - sum(b1^2*ipsi) -n*log(sigma^2) - pr*log(sigma.b^2) - 2*ldet - n*log(2*pi))/2
attr(l,"b") <- as.numeric(b1) ## return \hat beta and \hat b
-l
}

y<-matrix(c(110,100,110,100,100,110,110,100,100),nrow=9,ncol=1)
x<-matrix(c(1,1,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,1,1,1,1),nrow=9,ncol=3)

z<-matrix(c(1,0,0,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,0,
0,0,0,0,0,1,1,0,0,
0,1,0,1,1,0,0,1,1),nrow=9,ncol=4)

X1 <- cbind(x,z)

crossprod(X1)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 2 0 0 1 0 0 1
## [2,] 0 3 0 0 1 0 2
## [3,] 0 0 4 0 0 2 2
## [4,] 1 0 0 1 0 0 0
## [5,] 0 1 0 0 1 0 0
## [6,] 0 0 2 0 0 2 0
## [7,] 1 2 2 0 0 0 5
n <- length(y); pr <- ncol(z); pf <- ncol(x)

theta<-c(log(sqrt(0.1)),log(1))

sigma.b <- exp(theta[1])
sigma <- exp(theta[2])

ipsi <- c(rep(0,pf),rep(1/sigma.b^2,pr))

llm(theta,x,z,y)
## [1] 92.34967
## attr(,"b")
## [1] 105.6386574 104.2757378 105.4584366 0.3964857 0.5203875 0.7569272
## [7] -1.6738004

crossprod function

http://www.endmemo.com/program/R/crossprod.php

R crossprod Function

crossprod() function returns matrix cross-product.

crossprod(x, y = NULL)
tcrossprod(x, y = NULL)

x: numeric matrix
y: numeric matrix, if y=NULL, y is the same as x
...

> x <- matrix(1:9,3,3)
> x

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

> crossprod(x)

     [,1] [,2] [,3]
[1,]   14   32   50
[2,]   32   77  122
[3,]   50  122  194

> tcrossprod(x)

     [,1] [,2] [,3]
[1,]   66   78   90
[2,]   78   93  108
[3,]   90  108  126

»	R PCH Symbols
»	R Color Names
»	R Regular Expression
»	R tapply Function
»	R String Functions
»	R Plot Function
»	R Builtin Datasets List
	Python Tutorials
	HTML Tutorials
	JavaScript Tutorials
	Statistics
	News, Events Worldwide
	Unit Conversions
	Top Visited Websites Directory