R Under development (unstable) (2025-01-27 r87654 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #    Compare mmglm0 and mmglm1
> #    Gaussian with identity log function
> #    R CMD BATCH --no-save mmglm0-mmglm1-gaussian1.R mmglm0-mmglm1-gaussian1.Rout.save
> 
> library(HiddenMarkov)
> 
> 
> delta <- c(0,1)
> 
> Pi <- matrix(c(0.8, 0.2,
+                0.3, 0.7),
+              byrow=TRUE, nrow=2)
> 
> beta <- matrix(c(0.1, -0.1,
+                  1.0,  5.0),
+                byrow=TRUE, nrow=2)
> 
> sd <- c(1, 2)
> 
> #--------------------------------------------------------
> #     Gaussian with identity link function
> #         using mmglm0
> 
> x0 <- mmglm0(NULL, Pi, delta, family="gaussian", link="log",
+              beta=beta, sigma=sd, msg=FALSE)
> 
> x0 <- simulate(x0, nsim=5000, seed=10)
> 
> x0 <- BaumWelch(x0, bwcontrol(prt=FALSE))
> 
> print(summary(x0))
$variable.names
[1] "x1" "y" 

$delta
[1] 0 1

$Pi
          [,1]      [,2]
[1,] 0.8038487 0.1961513
[2,] 0.2841708 0.7158292

$nonstat
[1] TRUE

$beta
          [,1]       [,2]
[1,] 0.1210856 -0.1015701
[2,] 0.9762175  5.0001044

$sigma
[1] 0.989749 2.014274

$family
[1] "gaussian"

$glmformula
y ~ x1
<environment: 0x000001dc46c95450>

$link
[1] "log"

$n
[1] 5000

> 
> #--------------------------------------------------------
> #    Now embed this data into a mmglm1 object
> 
> glmformula <- formula(y ~ x1)
> glmfamily <- gaussian(link="log")
> Xdesign <- model.matrix(glmformula, data=x0$x)
> 
> x1 <- mmglm1(x0$x$y, Pi, delta, glmfamily, beta, Xdesign, sigma=sd, msg=FALSE)
> 
> x1 <- BaumWelch(x1, bwcontrol(prt=FALSE))
> 
> print(summary(x1))
$delta
[1] 0 1

$Pi
          [,1]      [,2]
[1,] 0.8038487 0.1961513
[2,] 0.2841708 0.7158292

$nonstat
[1] TRUE

$beta
              State 1    State 2
(Intercept) 0.1210856 -0.1015701
x1          0.9762175  5.0001044

$sigma
[1] 0.989749 2.014274

$glmfamily

Family: gaussian 
Link function: log 


$n
[1] 5000

> 
> #--------------------------------------------------------
> #   Compare Models
> 
> if (abs(logLik(x0)-logLik(x1)) > 1e-06)
+     warning("WARNING: See tests/mmglm0-mmglm1-gaussian1.R, log-likelihoods are different")
> 
> if (any(Viterbi(x0)!=Viterbi(x1)))
+     warning("WARNING: See tests/mmglm0-mmglm1-gaussian1.R, Viterbi paths are different")
> 
> if (any(abs(residuals(x0)-residuals(x1)) > 1e-06))
+     warning("WARNING: See tests/mmglm0-mmglm1-gaussian1.R, residuals are different")
> 
> 
> 
> proc.time()
   user  system elapsed 
   1.45    0.14    1.54