Curvas reales, media y simuladas

```{r} library(iAR) library(polynom) library(dplyr) library(readr) library(simex) library(xtable) require(Rdrw)

setwd(“C:/Users/cataluna/OneDrive/Escritorio/usach/Tesis/AGN/AGN”)

AGN1=read.table(“AGN1.csv”,sep=“,”,header=T) AGN2=read.table(“AGN2.csv”,sep=“,”,header=T) AGN3=read.table(“AGN3.csv”,sep=“,”,header=T) AGN4=read.table(“AGN4.csv”,sep=“,”,header=T) AGN5=read.table(“AGN5.csv”,sep=“,”,header=T) AGN6=read.table(“AGN6.csv”,sep=“,”,header=T) AGN7=read.table(“AGN7.csv”,sep=“,”,header=T) AGNF=rbind(AGN1,AGN2,AGN3,AGN4,AGN5,AGN6,AGN7)

aggregate(mjd~fid,AGNF,range)

object=names(table(AGNF[,1])) l0=length(object) crossmatchedobjects=read.csv(“ts_v9.0.1_SMBH_ZTF_xmatch.csv”) object=data.frame(“oid”=object) join=merge(object, crossmatchedobjects, by = “oid”)

MSEMLEg=rep(0,l0) o2=iAR::utilities() o2<-gentime(o2, n=400, distribution = “expmixture”, lambda1 = 55, lambda2 = 3,p1 = 0.15, p2 = 0.85)

newdata=NULL simdata=NULL

for(i in 1:l0) { obj=join$oid[i] class=join$survey_class_mapped[i] pos=which(AGNF[,1]==obj) lc=AGNF[pos,] lcg=na.omit(lc[which(lc[,3]==1),c(2,4,5)]) lcr=na.omit(lc[which(lc[,3]==2),c(2,4,5)]) lcg=lcg[order(lcg[,1]),] pos=which(diff(lcg[,1])==0)

if(length(pos)>0) { for(k in pos) { lcg[k,]=apply(lcg[c(k,k+1),],2,mean) } lcg=lcg[-c(pos+1),] }

p=which(lcg[,2]>=30)
if(length(p)>0) lcg=lcg[-p,]

if(dim(lcg)[1]<10) next

p=which(lcr[,2]>=30) if(length(p)>0) lcr=lcr[-p,]

lcr=lcr[order(lcr[,1]),] pos=which(diff(lcr[,1])==0)

if(length(pos)>0) { for(k in pos) { lcr[k,]=apply(lcr[c(k,k+1),],2,mean) } lcr=lcr[-c(pos+1),] }

if(dim(lcr)[1]<10) next

media=mean(lcg[,2]) primervalor=c(58242.22,media,0.01) ultimovalor=c(60316.16,media,0.01) lcg=rbind(primervalor,lcg,ultimovalor)

lcg <- lcg[order(lcg[,1]), , drop = FALSE]
lcg <- lcg[!duplicated(lcg[,1]), , drop = FALSE]
d <- diff(lcg[,1])

if (any(d <= 0)) {
eps <- .Machine$double.eps^0.5 for (j in which(d <= 0)) lcg[j+1,1] <- lcg[j,1] + eps }

fid=1 new=data.frame(obj,lcg[,1],fid,lcg[,2],lcg[,3]) names(new)=names(AGNF) newdata=rbind(newdata,new)

n=dim(lcg)[1] x <- iAR(family = “norm”, times = lcg[,1], series = lcg[,2], series_esd = lcg[,3]) x = loglik(x)

st <- sort(unique(as.numeric(o2@times)))
simiar <- iAR(family = “norm”, times = st, coef = x@coef) simiar <- iAR::sim(simiar) newsim <- data.frame(obj, st, fid, simiar@series + media, 0)

names(newsim)=names(AGNF) simdata=rbind(simdata,newsim) }

library(ggplot2) library(ggridges) library(patchwork)

purple_palette <- c( “AGN” = “#240046”,
“Blazar” = “#5A189A”,
“QSO” = “#E0AAFF”
)

panel_theme <- theme_minimal(base_size = 13) + theme( panel.border = element_rect(color = “black”, fill = NA, linewidth = 0.9), plot.margin = margin(6, 6, 6, 6), panel.grid.minor = element_blank() )

g1 <- ggplot( phi_no_out, aes(x = survey_class_mapped, y = phi, fill = survey_class_mapped) ) + geom_boxplot( color = “black”, alpha = 0.70, outlier.shape = NA, linewidth = 0.4, width = 0.6 ) + scale_fill_manual(values = purple_palette) + scale_y_continuous( limits = c(0.99, 1.0001), breaks = c(0.99, 0.995, 1.0), expand = expansion(mult = c(0.02, 0.08)) ) + labs( x = “Clase de objeto”, y = expression(phi) ) + panel_theme + theme( axis.text.x = element_text(size = 14, face = “bold”), axis.title.x = element_text(margin = margin(t = 6)), axis.title.y = element_text(margin = margin(r = 6)), legend.position = “none”
)

phi_plot <- phi_no_out phi_plot$survey_class_mapped <- reorder( phi_plot$survey_class_mapped, phi_plot$phi, FUN = median )

g2 <- ggplot( phi_plot, aes(x = phi, y = survey_class_mapped, fill = survey_class_mapped) ) + geom_density_ridges( color = “#2D004D”, alpha = 0.75, scale = 1.15, rel_min_height = 0.01, linewidth = 0.6 ) + scale_fill_manual(values = purple_palette) + labs( x = expression(phi), y = “Clase de objeto” ) + panel_theme + theme( axis.text.y = element_text(size = 14, face = “bold”), axis.title.x = element_text(margin = margin(t = 6)), axis.title.y = element_text(margin = margin(r = 6)), legend.position = “none”
)

final_plot <- (g1 | g2) + plot_layout(widths = c(1, 1)) + plot_annotation( title = “Distribución del coeficiente iAR (φ) por clase”, theme = theme( plot.title = element_text(hjust = 0.5, face = “bold”, size = 16), plot.margin = margin(6, 6, 6, 6) ) )

final_plot

library(dplyr) library(iAR) library(ggplot2)

cat(“Filas newdata:”, nrow(newdata), “”) cat(“Nombres newdata:”); print(names(newdata))

time_col <- if (“mjd” %in% names(newdata)) “mjd” else names(newdata)[2] mag_col <- if (“magpsf_corr” %in% names(newdata)) “magpsf_corr” else names(newdata)[4] esd_col <- if (“sigmapsf_corr” %in% names(newdata)) “sigmapsf_corr” else names(newdata)[5]

newdata[[time_col]] <- as.numeric(as.character(newdata[[time_col]])) newdata[[mag_col]] <- as.numeric(as.character(newdata[[mag_col]])) newdata[[esd_col]] <- as.numeric(as.character(newdata[[esd_col]]))

if (“fid” %in% names(newdata)) newdata$fid <- as.numeric(as.character(newdata$fid))

cat(“Valores únicos de fid en newdata:”) print(unique(newdata$fid))

use_fid <- 1 if (!any(newdata$fid == 1, na.rm = TRUE)) { use_fid <- sort(unique(na.omit(newdata$fid)))[1] cat(“OJO: no existe fid==1. Usaré fid =”, use_fid, “”) }

get_phi_from_fit <- function(fit) { cc <- fit@coef nms <- names(cc) if (!is.null(nms) && “phi” %in% nms) return(as.numeric(cc[“phi”])) if (!is.null(nms)) { idx <- which(grepl(“phi|ar”, nms, ignore.case = TRUE)) if (length(idx) > 0) return(as.numeric(cc[idx[1]])) } if (length(cc) > 0) return(as.numeric(cc[1])) NA_real_ }

calc_phi_one <- function(df) { t <- df[[time_col]] y <- df[[mag_col]] esd <- df[[esd_col]]

ok <- is.finite(t) & is.finite(y) & is.finite(esd) t <- t[ok]; y <- y[ok]; esd <- esd[ok] esd[!is.finite(esd) | esd <= 0] <- 0.01

o <- order(t) t <- t[o]; y <- y[o]; esd <- esd[o] dup <- duplicated(t) if (any(dup)) { t <- t[!dup]; y <- y[!dup]; esd <- esd[!dup] }

if (length(t) < 10) return(NA_real_)

tryCatch({ fit <- iAR(family=“norm”, times=t, series=y, series_esd=esd) fit <- loglik(fit) get_phi_from_fit(fit) }, error = function(e) { NA_real_ }) }

dat_f <- newdata %>% filter(fid == use_fid)

cat(“Filas después de filtrar fid:”, nrow(dat_f), “”) cat(“OIDs distintos:”, length(unique(dat_f$oid)), “”)

phi_df <- dat_f %>% group_by(oid) %>% summarise(phi = calc_phi_one(cur_data()), .groups = “drop”) %>% filter(is.finite(phi))

cat(“Phi calculados:”, nrow(phi_df), “”)

phi_class <- phi_df %>% left_join(join %>% select(oid, survey_class_mapped), by=“oid”) %>% filter(!is.na(survey_class_mapped))

cat(“Con clase asignada:”, nrow(phi_class), “”) print(table(phi_class$survey_class_mapped))

phi_no_out <- phi_class %>% group_by(survey_class_mapped) %>% filter(n() >= 5) %>% mutate( q1 = quantile(phi, 0.25, na.rm=TRUE), q3 = quantile(phi, 0.75, na.rm=TRUE), iqr = q3 - q1, low = q1 - 1.5iqr, high = q3 + 1.5iqr ) %>% filter(phi >= low, phi <= high) %>% ungroup()

cat(“Post-outliers:”, nrow(phi_no_out), “”)

if (nrow(phi_no_out) == 0) { cat(“No hay datos para graficar tras IQR. Probemos graficar sin filtrar outliers.”) } else { ggplot(phi_no_out, aes(x = survey_class_mapped, y = phi)) + geom_boxplot(fill=“#7B2CBF”, color=“black”, alpha=0.35, outlier.shape=NA) + scale_y_continuous(limits=c(0.99, 1.0001), breaks=c(0.99, 0.995, 1.0)) + labs(title=“Distribución del Coeficiente iAR (phi) por Clase”, x=“Clase de Objeto”, y=“Coeficiente (phi)”) + theme_minimal(base_size=12) + theme( plot.title = element_text(hjust=0.5, face=“bold”), panel.grid.major = element_line(color=“grey85”), panel.grid.minor = element_line(color=“grey92”) ) }

library(dplyr) library(iAR) library(ggplot2) library(ggridges) library(patchwork)

cat(“Filas newdata:”, nrow(newdata), “”) cat(“Nombres newdata:”); print(names(newdata))

if (“fid” %in% names(newdata)) newdata$fid <- as.numeric(as.character(newdata$fid))

cat(“Valores únicos de fid en newdata:”) print(unique(newdata$fid))

use_fid <- 1 if (!any(newdata$fid == 1, na.rm = TRUE)) { use_fid <- sort(unique(na.omit(newdata$fid)))[1] cat(“OJO: no existe fid==1. Usaré fid =”, use_fid, “”) }

calc_phi_one <- function(df) { t <- df[[time_col]] y <- df[[mag_col]] esd <- df[[esd_col]]

ok <- is.finite(t) & is.finite(y) & is.finite(esd) t <- t[ok]; y <- y[ok]; esd <- esd[ok] esd[!is.finite(esd) | esd <= 0] <- 0.01

o <- order(t) t <- t[o]; y <- y[o]; esd <- esd[o] dup <- duplicated(t) if (any(dup)) { t <- t[!dup]; y <- y[!dup]; esd <- esd[!dup] }

if (length(t) < 10) return(NA_real_)

tryCatch({ fit <- iAR(family=“norm”, times=t, series=y, series_esd=esd) fit <- loglik(fit) get_phi_from_fit(fit) }, error = function(e) NA_real_) }

dat_f <- newdata %>% filter(fid == use_fid)

cat(“Filas después de filtrar fid:”, nrow(dat_f), “”) cat(“OIDs distintos:”, length(unique(dat_f$oid)), “”)

phi_df <- dat_f %>% group_by(oid) %>% summarise(phi = calc_phi_one(cur_data()), .groups = “drop”) %>% filter(is.finite(phi))

cat(“Phi calculados:”, nrow(phi_df), “”)

phi_class <- phi_df %>% left_join(join %>% select(oid, survey_class_mapped), by=“oid”) %>% filter(!is.na(survey_class_mapped))

cat(“Con clase asignada:”, nrow(phi_class), “”) print(table(phi_class$survey_class_mapped))

cat(“Post-outliers:”, nrow(phi_no_out), “”)

if (nrow(phi_no_out) == 0) { cat(“No hay datos para graficar tras IQR.”) } else {

g1 <- ggplot(phi_no_out, aes(x = survey_class_mapped, y = phi)) + geom_boxplot(fill=“#7B2CBF”, color=“black”, alpha=0.35, outlier.shape=NA, linewidth=0.4) + scale_y_continuous(limits=c(0.998, 1.0001), breaks=c(0.998, 0.999, 1.0)) + labs(title=“(a) Boxplot de 3C6 por clase”, x=“Clase de Objeto”, y=expression(phi)) + theme_minimal(base_size=12) + theme( plot.title = element_text(hjust=0.5, face=“bold”), axis.text.x = element_text(size = 14, face=“bold”), panel.grid.major = element_line(color=“grey85”), panel.grid.minor = element_line(color=“grey92”) )

phi_plot <- phi_no_out phi_plot$survey_class_mapped <- reorder(phi_plot$survey_class_mapped, phi_plot$phi, FUN = median)

g2 <- ggplot(phi_plot, aes(x = phi, y = survey_class_mapped)) + geom_density_ridges( fill = “#7B2CBF”, alpha = 0.60, color = “#2D004D”, scale = 1.2, rel_min_height = 0.01 ) + labs( title = “(b) Densidad de 3C6 por clase”, x = expression(phi), y = “Clase de objeto” ) + theme_minimal(base_size = 14) + theme( plot.title = element_text(hjust = 0.5, face = “bold”, size = 14), axis.text.y = element_text(size = 14, face = “bold”) # <- agranda clases )

(g1 | g2) + plot_layout(widths = c(1, 1.15)) }

library(ggplot2) library(ggridges) library(patchwork)

g1 <- ggplot( phi_no_out, aes(x = survey_class_mapped, y = phi, fill = survey_class_mapped) ) + geom_boxplot( color = “black”, alpha = 0.70, outlier.shape = NA, linewidth = 0.4, width = 0.6 ) + scale_fill_manual(values = purple_palette) + scale_y_continuous( limits = c(0.998, 1.0001), breaks = c(0.998, 0.999, 1.0), expand = expansion(mult = c(0.02, 0.08)) ) + labs( x = “Clase de objeto”, y = expression(phi) ) + panel_theme + theme( axis.text.x = element_text(size = 14, face = “bold”), axis.title.x = element_text(margin = margin(t = 6)), axis.title.y = element_text(margin = margin(r = 6)), legend.position = “none” )

phi_plot <- phi_no_out phi_plot$survey_class_mapped <- reorder( phi_plot$survey_class_mapped, phi_plot$phi, FUN = median )

final_plot

if (!exists(“newdata”)) stop(“No existe ‘newdata’. Ejecuta primero el loop.”)

primeros_obj <- unique(newdata$oid)[1:3]

data_plot <- newdata[newdata$oid %in% primeros_obj, ]

par(mfrow = c(3,1), mar = c(4,4,2,1))

for(obj in primeros_obj){

d <- data_plot[data_plot$oid == obj, ] d <- d[order(d$mjd), ]

plot( d$mjd, d$magpsf_corr, type = “l”, main = paste(“Objeto:”, obj), xlab = “Tiempo (MJD)”, ylab = “Magnitud”, col = “black”, lwd = 1.6 )

points( d$mjd, d$magpsf_corr, pch = 20, col = “black”, cex = 0.8 ) }

Volver a layout normal

par(mfrow = c(1,1))

if (!exists(“newdata”)) stop(“No existe ‘newdata’. Ejecuta primero el loop.”)

primeros_obj <- unique(AGNF$oid)[1:3] data_plot <- newdata[AGNF$oid %in% primeros_obj, ]

par(mfrow = c(3,1), mar = c(4,4,2,1), oma = c(0,0,3,0))

for(obj in primeros_obj){

d <- data_plot[data_plot$oid == obj, ] d <- d[order(d$mjd), ]

plot( d$mjd, d$magpsf_corr, type = “l”, main = paste(“Objeto:”, obj), xlab = “Tiempo (MJD)”, ylab = “Magnitud”, col = “black”, lwd = 1.7 )

points( d$mjd, d$magpsf_corr, pch = 20, col = “black”, cex = 0.8 ) }

Título general

mtext(“Primeras tres curvas de luz reales (banda g)”, outer = TRUE, cex = 1, font = 2)

par(mfrow = c(1,1))

Mostrar primeras 3 curvas ORIGINALES (sin corrección artificial)

primeros_obj <- unique(AGNF$oid)[1:3]

par(mfrow = c(3,1), mar = c(4,4,2,1), oma = c(0,0,3,0))

for(obj in primeros_obj){

lc <- AGNF[AGNF$oid == obj, ]

lc <- lc[lc$fid == 1, c(“mjd”,“magpsf_corr”,“sigmapsf_corr_ext”)] lc <- na.omit(lc)

lc <- lc[order(lc$mjd), ]

lc <- lc[!duplicated(lc$mjd), ]

if(nrow(lc) < 2) next

plot( lc$mjd, lc$magpsf_corr, type = “l”, main = paste(“Objeto:”, obj), xlab = “Tiempo (MJD)”, ylab = “Magnitud”, col = “black”, lwd = 1.7 )

points( lc$mjd, lc$magpsf_corr, pch = 20, col = “black”, cex = 0.8 ) }

mtext(“Primeras tres curvas de luz (banda g)”, outer = TRUE, cex = 1.3, font = 2)

par(mfrow = c(1,1))