diff --git a/NAMESPACE b/NAMESPACE
index 8cd36c9..de50687 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(coef,ggwrm)
 S3method(coef,gtdrm)
 S3method(coef,gwlcrm)
 S3method(coef,gwrm)
@@ -9,6 +10,7 @@ S3method(fitted,gwlcrm)
 S3method(fitted,gwrm)
 S3method(fitted,gwrmultiscalem)
 S3method(glyph.plot,gwpcam)
+S3method(plot,ggwrm)
 S3method(plot,gtdrm)
 S3method(plot,gwcorrm)
 S3method(plot,gwlcrm)
@@ -16,6 +18,7 @@ S3method(plot,gwrm)
 S3method(plot,gwrmultiscalem)
 S3method(plot,modelselcritl)
 S3method(predict,gwrm)
+S3method(print,ggwrm)
 S3method(print,gtdrm)
 S3method(print,gwavgm)
 S3method(print,gwcorrm)
@@ -38,6 +41,7 @@ export(gwcorr_config)
 export(gwcorrelation)
 export(gwpca)
 export(gwr_basic)
+export(gwr_generalized)
 export(gwr_lcr)
 export(gwr_multiscale)
 export(mgwr_config)
diff --git a/R/RcppExports.R b/R/RcppExports.R
index d783a29..61ea322 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -25,6 +25,14 @@ gwr_basic_predict <- function(pcoords, x, y, coords, bw, adaptive, kernel, longl
     .Call(`_GWmodel3_gwr_basic_predict`, pcoords, x, y, coords, bw, adaptive, kernel, longlat, p, theta, intercept, parallel_type, parallel_arg, verbose)
 }
 
+gwr_generalized_fit <- function(x, y, coords, family, bw, adaptive, kernel, longlat, p, theta, hatmatrix, intercept, optim_bw, optim_bw_criterion, parallel_type, parallel_arg) {
+    .Call(`_GWmodel3_gwr_generalized_fit`, x, y, coords, family, bw, adaptive, kernel, longlat, p, theta, hatmatrix, intercept, optim_bw, optim_bw_criterion, parallel_type, parallel_arg)
+}
+
+gwr_generalized_predict <- function(pcoords, x, y, coords, family, bw, adaptive, kernel, longlat, p, theta, hatmatrix, intercept, parallel_type, parallel_arg) {
+    .Call(`_GWmodel3_gwr_generalized_predict`, pcoords, x, y, coords, family, bw, adaptive, kernel, longlat, p, theta, hatmatrix, intercept, parallel_type, parallel_arg)
+}
+
 gwr_lcr_fit <- function(x, y, coords, bw, adaptive, kernel, longlat, p, theta, lambda, lambda_adjust, cn_thresh, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw) {
     .Call(`_GWmodel3_gwr_lcr_fit`, x, y, coords, bw, adaptive, kernel, longlat, p, theta, lambda, lambda_adjust, cn_thresh, intercept, hatmatrix, parallel_type, parallel_arg, optim_bw)
 }
diff --git a/R/gwr_generalized.R b/R/gwr_generalized.R
new file mode 100644
index 0000000..1113eac
--- /dev/null
+++ b/R/gwr_generalized.R
@@ -0,0 +1,371 @@
+#' Calibrate a basic GGWR model
+#'
+#' @param formula Regresison model.
+#' @param data A `sf` objects.
+#' @param bw Either a value to set the size of bandwidth,
+#'  or one of the following characters to set the criterion for
+#'  bandwidth auto-optimization process.
+#'  - `AIC`
+#'  - `CV`
+#'  Note that if `NA` or other non-numeric value is setted,
+#'  this parameter will be reset to `Inf`.
+#' @param adaptive Whether the bandwidth value is adaptive or not.
+#' @param kernel Kernel function used.
+#' @param longlat Whether the coordinates
+#' @param p Power of the Minkowski distance,
+#'  default to 2, i.e., Euclidean distance.
+#' @param theta Angle in radian to roate the coordinate system, default to 0.
+#' @param hatmatrix If TRUE, great circle will be caculated.
+#' @param parallel_method Parallel method.
+#' @param parallel_arg Parallel method argument.
+#' @param verbose Whether to print additional information.
+#'
+#' @return A `ggwrm` object.
+#' 
+#' @details
+#' ## Parallelization
+#' 
+#' One parallel methods are provided to speed up basic GGWR algorithm:
+#' 
+#' - Multithreading (`omp`)
+#' 
+#' See the vignettes about parallelization to learn more about this topic.
+#'
+#' @examples
+#' data(LondonHP)
+#'
+#' # Basic usage
+#' gwr_generalized(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, "poisson", 64, TRUE)
+#'
+#' # Bandwidth Optimization
+#' m <- gwr_generalized(PURCHASE ~ FLOORSZ + UNEMPLOY + PROF, LondonHP, 'AIC', TRUE)
+#' m
+#' 
+#' @seealso `browseVignettes("")`
+#'
+#' @importFrom stats na.action model.frame model.extract model.matrix terms
+#' @export
+gwr_generalized <- function(
+    formula,
+    data,
+    family = c("poisson", "binomial"),
+    bw = NA,
+    adaptive = FALSE,
+    kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"),
+    longlat = FALSE,
+    p = 2.0,
+    theta = 0.0,
+    hatmatrix = TRUE,
+    parallel_method = c("no", "omp"),
+    parallel_arg = c(0)
+) {
+    ### Check args
+    if(missing(family)){
+        family = "poisson"
+    }
+    family = match.arg(family)
+    kernel = match.arg(kernel)
+    parallel_method = match.arg(parallel_method)
+    attr(data, "na.action") <- getOption("na.action")
+
+    ### Extract coords
+    data <- do.call(na.action(data), args = list(data))
+    coords <- as.matrix(sf::st_coordinates(sf::st_centroid(data)))
+    if (is.null(coords) || nrow(coords) != nrow(data))
+        stop("Missing coordinates.")
+
+    ### Extract variables
+    mc <- match.call(expand.dots = FALSE)
+    mf <- model.frame(formula = formula(formula), data = sf::st_drop_geometry(data))
+    mt <- attr(mf, "terms")
+    y <- model.extract(mf, "response")
+    x <- model.matrix(mt, mf)
+    dep_var <- as.character(attr(terms(mf), "variables")[[2]])
+    has_intercept <- attr(terms(mf), "intercept") == 1
+    indep_vars <- colnames(x)
+    indep_vars[which(indep_vars == "(Intercept)")] <- "Intercept"
+    colnames(x) <- indep_vars
+    if (has_intercept && indep_vars[1] != "Intercept") {
+        stop("Please put Intercept to the first column.")
+    }
+
+    ### Check whether bandwidth is valid.
+    if (missing(bw)) {
+        optim_bw <- TRUE
+        optim_bw_criterion <- "AIC"
+        bw <- Inf
+    } else if (is.numeric(bw) || is.integer(bw)) {
+        optim_bw <- FALSE
+        optim_bw_criterion <- "AIC"
+    } else {
+        optim_bw <- TRUE
+        optim_bw_criterion <-
+            ifelse(is.character(bw), match.arg(bw, c("CV", "AIC")), "AIC")
+        bw <- Inf
+    }
+
+    ### Call solver
+    c_result <- tryCatch(gwr_generalized_fit(
+        x, y, coords, 
+        enum(family, c("poisson", "binomial")),
+        bw, adaptive, enum(kernel), longlat, p, theta,
+        hatmatrix, has_intercept,
+        optim_bw, enum(optim_bw_criterion, c("AIC", "CV")),
+        enum_list(parallel_method, parallel_types), parallel_arg
+    ), error = function (e) {
+        stop("Error:", conditionMessage(e))
+    })
+    if (optim_bw)
+        bw <- c_result$bandwidth
+    betas <- c_result$betas
+    diagnostic <- c_result$diagnostic
+
+    yhat <- c_result$yhat
+    resi <- y - yhat
+    # n_dp <- nrow(coords)
+    # rss_gw <- sum(resi * resi)
+    # sigma <- rss_gw / (n_dp - 2 * shat_trace[1] + shat_trace[2])
+    # betas_se <- sqrt(sigma * betas_se)
+    # betas_tv <- betas / betas_se
+
+    ### Create result Layer
+    colnames(betas) <- indep_vars
+    # colnames(betas_se) <- paste(indep_vars, "SE", sep = ".")
+    # colnames(betas_tv) <- paste(indep_vars, "TV", sep = ".")
+    sdf_data <- data.frame(
+        betas,
+        yhat = yhat,
+        residual = resi
+        # betas_se,
+        # betas_tv
+    )
+    sdf_data$geometry <- sf::st_geometry(data)
+    sdf <- sf::st_sf(sdf_data)
+
+    ### Return result
+    ggwrm <- list(
+        SDF = sdf,
+        diagnostic = diagnostic,
+        args = list(
+            x = x,
+            y = y,
+            coords = coords,
+            family = family,
+            bw = bw,
+            adaptive = adaptive,
+            kernel = kernel,
+            longlat = longlat,
+            p = p,
+            theta = theta,
+            hatmatrix = hatmatrix,
+            has_intercept = has_intercept,
+            optim_bw = optim_bw,
+            optim_bw_criterion = optim_bw_criterion,
+            parallel_method = parallel_method,
+            parallel_arg = parallel_arg
+        ),
+        call = mc,
+        indep_vars = indep_vars,
+        dep_var = dep_var
+    )
+    class(ggwrm) <- "ggwrm"
+    ggwrm
+}
+
+#' Print description of a `ggwrm` object
+#'
+#' @param x An `ggwrm` object returned by [gwr_generalized()].
+#' @param decimal_fmt The format string passing to [base::sprintf()].
+#' @inheritDotParams print_table_md
+#' 
+#' @method print ggwrm
+#' @importFrom stats coef fivenum
+#' @rdname print
+#' @export
+print.ggwrm <- function(x, decimal_fmt = "%.3f", ...) {
+    if (!inherits(x, "ggwrm")) {
+        stop("It's not a ggwrm object.")
+    }
+
+    ### Basic Information
+    cat("Generalized Geographically Weighted Regression Model", fill = T)
+    cat("========================================", fill = T)
+    cat("  Formula:", deparse(x$call$formula), fill = T)
+    cat("     Data:", deparse(x$call$data), fill = T)
+    cat("   Family:", x$args$family, fill = T)
+    cat("   Kernel:", x$args$kernel, fill = T)
+    cat("Bandwidth:", x$args$bw,
+        ifelse(x$args$adaptive, "(Nearest Neighbours)", "(Meters)"),
+        ifelse(x$args$optim_bw, paste0(
+            "(Optimized accroding to ",
+            x$args$optim_bw_criterion,
+            ")"
+        ), ""), fill = T)
+    cat("\n", fill = T)
+
+    cat("Summary of Coefficient Estimates", fill = T)
+    cat("--------------------------------", fill = T)
+    betas <- coef(x)
+    beta_fivenum <- t(apply(betas, 2, fivenum))
+    colnames(beta_fivenum) <- c("Min.", "1st Qu.", "Median", "3rd Qu.", "Max.")
+    rownames(beta_fivenum) <- colnames(betas)
+    beta_str <- rbind(
+        c("Coefficient", colnames(beta_fivenum)),
+        cbind(rownames(beta_fivenum), matrix2char(beta_fivenum, decimal_fmt))
+    )
+    print_table_md(beta_str, ...)
+    cat("\n", fill = T)
+
+    cat("Diagnostic Information", fill = T)
+    cat("----------------------", fill = T)
+    cat("  RSS:", x$diagnostic$RSS, fill = T)
+    cat("   R2:", x$diagnostic$RSquare, fill = T)
+    cat("  AIC:", x$diagnostic$AIC, fill = T)
+    cat(" AICc:", x$diagnostic$AICc, fill = T)
+    cat("\n", fill = T)
+}
+
+#' @describeIn gwr_generalized Plot the result of basic GGWR model.
+#'
+#' @param x A "ggwrm" object.
+#' @param y Ignored.
+#' @param columns Column names to plot.
+#'  If it is missing or non-character value, all coefficient columns are plottd.
+#' @param \dots Additional arguments passing to [sf::plot()].
+#' @method plot ggwrm
+#'
+#' @examples
+#' plot(m)
+#'
+#' @export
+plot.ggwrm <- function(x, y, ..., columns) {
+    if (!inherits(x, "ggwrm")) {
+        stop("It's not a ggwrm object.")
+    }
+
+    sdf <- sf::st_as_sf(x$SDF)
+    sdf_colnames <- names(sf::st_drop_geometry(x$SDF))
+    if (!missing(columns) && is.character(columns)) {
+        valid_columns <- intersect(columns, sdf_colnames)
+        if (length(valid_columns) > 0) {
+            sdf <- sdf[valid_columns]
+        }
+    } else { ### Select coefficient columns.
+        sdf <- sdf[x$indep_vars]
+    }
+    plot(sdf, ...)
+}
+
+#' @describeIn gwr_generalized Get coefficients of a Generalized GWR model.
+#'
+#' @param object A "ggwrm" object.
+#' @param \dots Additional arguments passing to [coef()].
+#' 
+#' @examples
+#' coef(m)
+#' 
+#' @method coef ggwrm
+#' @export
+coef.ggwrm <- function(object, ...) {
+    if (!inherits(object, "ggwrm")) {
+        stop("It's not a ggwrm object.")
+    }
+    sf::st_drop_geometry(object$SDF[object$indep_vars])
+}
+
+# ### 测试问题可能在内核库的predict函数中，应该先设置一下mHasRegressionData，否则导致现在出现矩阵列数计算的错误。
+# ### 在内核库中测试通过predict后，再完善这部分内容。
+# #' @describeIn gwr_generalized Predict on new locations.
+# #'
+# #' @param object A "ggwrm" object.
+# #' @param regression_points Data of new locations.
+# #' @param \dots Additional arguments.
+# #' @param verbose Whether to print additional message.
+# #'
+# #' @method predict ggwrm
+# #'
+# #' @examples
+# #' ### predict(m, LondonHP)
+# #'
+# #' @export
+# predict.ggwrm <- function(object, regression_points, verbose = FALSE, ...) {
+#     if (!inherits(object, "ggwrm")) {
+#         stop("It's not a ggwrm object.")
+#     }
+#     if (!inherits(regression_points, c("sf", "sfc", "matrix", "data.frame"))) {
+#         stop("The type of regression_points is not supported.")
+#     }
+
+#     has_intercept <- object$args$has_intercept
+#     ### Extract coordinates (and data, if possible)
+#     pcoords <- NULL
+#     py <- NULL
+#     px <- NULL
+#     if (inherits(regression_points, c("sf", "sfc"))) {
+#         pcoords <- sf::st_coordinates(sf::st_centroid(regression_points))
+#         ### Extract X
+#         indep_vars <- object$indep_vars
+#         if (has_intercept) indep_vars <- indep_vars[-1]
+#         if (all(indep_vars %in% names(regression_points))) {
+#             px <- as.matrix(sf::st_drop_geometry(regression_points[indep_vars]))
+#         }
+#         if (has_intercept) {
+#             px <- cbind(1, px)
+#         }
+#         ### Extract Y
+#         dep_var <- object$dep_var
+#         if (dep_var %in% names(regression_points)) {
+#             py <- regression_points[[dep_var]]
+#         }
+#     } else if(inherits(regression_points, c("matrix", "data.frame"))) {
+#         if (length(dim(regression_points)) == 2 && ncol(regression_points) == 2) {
+#             pcoords <- as.matrix(regression_points)
+#             if (!is.numeric(pcoords)) {
+#                 stop("Cannot convert regression_points to numeric matrix.")
+#             }
+#         } else {
+#             stop("Only matrix or data.frame of 2 columns is supported.")
+#         }
+#     } else {
+#         stop("Cannot extract coordinates of regression_points")
+#     }
+#     pcoords <- as.matrix(pcoords)
+
+#     ### refer to libgwmodel code GWRGenetalized.cpp:175
+#     if(object$args$hatmatrix) warning("hatmatrix is not need in prediction")
+#     has_hatmatrix <- FALSE
+
+#     ### Predict coefficients
+#     c_betas <- tryCatch(with(object$args, gwr_generalized_predict(
+#         pcoords, x, y, coords,
+#         enum(family, c("poisson", "binomial")),
+#         bw, adaptive, enum(kernel, kernel_enums),
+#         longlat, p, theta,
+#         has_hatmatrix, has_intercept,
+#         enum_list(parallel_method, parallel_types), parallel_arg
+#     )), error = function(e) {
+#         stop("Error:", conditionMessage(e))
+#     })
+
+#     result <- as.data.frame(c_betas)
+#     colnames(result) <- object$indep_vars
+
+#     ### If px exists, calculate yhat
+#     if (!is.null(px)) {
+#         yhat <- rowSums(px * c_betas)
+#         result <- cbind(result, "yhat" = yhat)
+
+#         if (!is.null(py)) {
+#             resi <- py - yhat
+#             result <- cbind(result, "residual" = resi)
+#         }
+#     }
+
+#     ### If regression_points is a sf object, return sf
+#     if (inherits(regression_points, c("sf", "sfc"))) {
+#         result$geometry <- sf::st_geometry(regression_points)
+#         result <- sf::st_as_sf(result)
+#     }
+
+#     result
+# }
\ No newline at end of file
diff --git a/man/gwr_generalized.Rd b/man/gwr_generalized.Rd
new file mode 100644
index 0000000..b930d2f
--- /dev/null
+++ b/man/gwr_generalized.Rd
@@ -0,0 +1,114 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gwr_generalized.R
+\name{gwr_generalized}
+\alias{gwr_generalized}
+\alias{plot.ggwrm}
+\alias{coef.ggwrm}
+\title{Calibrate a basic GGWR model}
+\usage{
+gwr_generalized(
+  formula,
+  data,
+  family = c("poisson", "binomial"),
+  bw = NA,
+  adaptive = FALSE,
+  kernel = c("gaussian", "exp", "bisquare", "tricube", "boxcar"),
+  longlat = FALSE,
+  p = 2,
+  theta = 0,
+  hatmatrix = TRUE,
+  parallel_method = c("no", "omp"),
+  parallel_arg = c(0)
+)
+
+\method{plot}{ggwrm}(x, y, ..., columns)
+
+\method{coef}{ggwrm}(object, ...)
+}
+\arguments{
+\item{formula}{Regresison model.}
+
+\item{data}{A \code{sf} objects.}
+
+\item{bw}{Either a value to set the size of bandwidth,
+or one of the following characters to set the criterion for
+bandwidth auto-optimization process.
+\itemize{
+\item \code{AIC}
+\item \code{CV}
+Note that if \code{NA} or other non-numeric value is setted,
+this parameter will be reset to \code{Inf}.
+}}
+
+\item{adaptive}{Whether the bandwidth value is adaptive or not.}
+
+\item{kernel}{Kernel function used.}
+
+\item{longlat}{Whether the coordinates}
+
+\item{p}{Power of the Minkowski distance,
+default to 2, i.e., Euclidean distance.}
+
+\item{theta}{Angle in radian to roate the coordinate system, default to 0.}
+
+\item{hatmatrix}{If TRUE, great circle will be caculated.}
+
+\item{parallel_method}{Parallel method.}
+
+\item{parallel_arg}{Parallel method argument.}
+
+\item{x}{A "ggwrm" object.}
+
+\item{y}{Ignored.}
+
+\item{\dots}{Additional arguments passing to \code{\link[=coef]{coef()}}.}
+
+\item{columns}{Column names to plot.
+If it is missing or non-character value, all coefficient columns are plottd.}
+
+\item{object}{A "ggwrm" object.}
+
+\item{verbose}{Whether to print additional information.}
+}
+\value{
+A \code{ggwrm} object.
+}
+\description{
+Calibrate a basic GGWR model
+}
+\details{
+\subsection{Parallelization}{
+
+One parallel methods are provided to speed up basic GGWR algorithm:
+\itemize{
+\item Multithreading (\code{omp})
+}
+
+See the vignettes about parallelization to learn more about this topic.
+}
+}
+\section{Functions}{
+\itemize{
+\item \code{plot(ggwrm)}: Plot the result of basic GGWR model.
+
+\item \code{coef(ggwrm)}: Get coefficients of a Generalized GWR model.
+
+}}
+\examples{
+data(LondonHP)
+
+# Basic usage
+gwr_generalized(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, "poisson", 64, TRUE)
+
+# Bandwidth Optimization
+m <- gwr_generalized(PURCHASE ~ FLOORSZ + UNEMPLOY + PROF, LondonHP, 'AIC', TRUE)
+m
+
+plot(m)
+
+coef(m)
+
+}
+\seealso{
+\code{browseVignettes("")}
+}
diff --git a/man/print.Rd b/man/print.Rd
index d55484e..95a0c3c 100644
--- a/man/print.Rd
+++ b/man/print.Rd
@@ -1,12 +1,14 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/gtdr.R, R/gwaverage.R, R/gwcorrelation.R,
-%   R/gwpca.R, R/gwr_basic.R, R/gwr_localcollinearity.R, R/gwr_multiscale.R
+%   R/gwpca.R, R/gwr_basic.R, R/gwr_generalized.R, R/gwr_localcollinearity.R,
+%   R/gwr_multiscale.R
 \name{print.gtdrm}
 \alias{print.gtdrm}
 \alias{print.gwavgm}
 \alias{print.gwcorrm}
 \alias{print.gwpcam}
 \alias{print.gwrm}
+\alias{print.ggwrm}
 \alias{print.gwlcrm}
 \alias{print.gwrmultiscalem}
 \title{Print description of a \code{gtdrm} object}
@@ -21,6 +23,8 @@
 
 \method{print}{gwrm}(x, decimal_fmt = "\%.3f", ...)
 
+\method{print}{ggwrm}(x, decimal_fmt = "\%.3f", ...)
+
 \method{print}{gwlcrm}(x, decimal_fmt = "\%.3f", ...)
 
 \method{print}{gwrmultiscalem}(x, decimal_fmt = "\%.3f", ...)
@@ -31,7 +35,7 @@
 \item{decimal_fmt}{The format string passing to \code{\link[base:sprintf]{base::sprintf()}}.}
 
 \item{...}{
-  Arguments passed on to \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}
+  Arguments passed on to \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}, \code{\link[=print_table_md]{print_table_md}}
   \describe{
     \item{\code{col.sep}}{Column seperator. Default to \code{""}.}
     \item{\code{header.sep}}{Header seperator. Default to \code{"-"}.}
@@ -54,6 +58,8 @@ Print description of a \code{gwpcam} object
 
 Print description of a \code{gwrm} object
 
+Print description of a \code{ggwrm} object
+
 Print description of a \code{gwlcrm} object
 
 Print description of a \code{gwrmultiscalem} object
diff --git a/src/Makevars.in b/src/Makevars.in
index 75f73ed..9201bef 100644
--- a/src/Makevars.in
+++ b/src/Makevars.in
@@ -41,9 +41,13 @@ OBJECTS_LIBGWMODEL = \
 	libgwmodel/src/gwmodelpp/spatialweight/SpatialWeight.o \
 	libgwmodel/src/gwmodelpp/spatialweight/Weight.o \
 	libgwmodel/src/gwmodelpp/BandwidthSelector.o \
+	libgwmodel/src/gwmodelpp/GeneralizedLinearModel.o \
+	libgwmodel/src/gwmodelpp/PoissonModel.o \
+	libgwmodel/src/gwmodelpp/BinomialModel.o \
 	libgwmodel/src/gwmodelpp/GTDR.o \
 	libgwmodel/src/gwmodelpp/GWRBase.o \
 	libgwmodel/src/gwmodelpp/GWRBasic.o \
+	libgwmodel/src/gwmodelpp/GWRGeneralized.o \
 	libgwmodel/src/gwmodelpp/GWRMultiscale.o \
 	libgwmodel/src/gwmodelpp/GWAverage.o \
 	libgwmodel/src/gwmodelpp/GWCorrelation.o \
@@ -64,6 +68,7 @@ OBJECTS_TELEGRAM = \
 OBJECTS_GWMODEL = \
 	utils.o \
 	gwr_basic.o \
+	gwr_generalized.o \
 	gwr_multiscale.o \
 	gwr_localcollinearity.o \
 	gtdr.o \
diff --git a/src/Makevars.win b/src/Makevars.win
index 7f4fcce..f657e3b 100644
--- a/src/Makevars.win
+++ b/src/Makevars.win
@@ -43,9 +43,13 @@ OBJECTS_LIBGWMODEL = \
 	libgwmodel/src/gwmodelpp/spatialweight/SpatialWeight.o \
 	libgwmodel/src/gwmodelpp/spatialweight/Weight.o \
 	libgwmodel/src/gwmodelpp/BandwidthSelector.o \
+	libgwmodel/src/gwmodelpp/GeneralizedLinearModel.o \
+	libgwmodel/src/gwmodelpp/PoissonModel.o \
+	libgwmodel/src/gwmodelpp/BinomialModel.o \
 	libgwmodel/src/gwmodelpp/GTDR.o \
 	libgwmodel/src/gwmodelpp/GWRBase.o \
 	libgwmodel/src/gwmodelpp/GWRBasic.o \
+	libgwmodel/src/gwmodelpp/GWRGeneralized.o \
 	libgwmodel/src/gwmodelpp/GWRMultiscale.o \
 	libgwmodel/src/gwmodelpp/GWAverage.o \
 	libgwmodel/src/gwmodelpp/GWCorrelation.o \
@@ -66,6 +70,7 @@ OBJECTS_TELEGRAM = \
 OBJECTS_GWMODEL = \
 	utils.o \
 	gwr_basic.o \
+	gwr_generalized.o \
 	gwr_multiscale.o \
 	gwr_localcollinearity.o \
 	gtdr.o \
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 02a2a2a..91f4dc0 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -161,6 +161,57 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// gwr_generalized_fit
+List gwr_generalized_fit(const arma::mat& x, const arma::vec& y, const arma::mat& coords, size_t family, double bw, bool adaptive, size_t kernel, bool longlat, double p, double theta, bool hatmatrix, bool intercept, bool optim_bw, size_t optim_bw_criterion, size_t parallel_type, const IntegerVector& parallel_arg);
+RcppExport SEXP _GWmodel3_gwr_generalized_fit(SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP familySEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP hatmatrixSEXP, SEXP interceptSEXP, SEXP optim_bwSEXP, SEXP optim_bw_criterionSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type coords(coordsSEXP);
+    Rcpp::traits::input_parameter< size_t >::type family(familySEXP);
+    Rcpp::traits::input_parameter< double >::type bw(bwSEXP);
+    Rcpp::traits::input_parameter< bool >::type adaptive(adaptiveSEXP);
+    Rcpp::traits::input_parameter< size_t >::type kernel(kernelSEXP);
+    Rcpp::traits::input_parameter< bool >::type longlat(longlatSEXP);
+    Rcpp::traits::input_parameter< double >::type p(pSEXP);
+    Rcpp::traits::input_parameter< double >::type theta(thetaSEXP);
+    Rcpp::traits::input_parameter< bool >::type hatmatrix(hatmatrixSEXP);
+    Rcpp::traits::input_parameter< bool >::type intercept(interceptSEXP);
+    Rcpp::traits::input_parameter< bool >::type optim_bw(optim_bwSEXP);
+    Rcpp::traits::input_parameter< size_t >::type optim_bw_criterion(optim_bw_criterionSEXP);
+    Rcpp::traits::input_parameter< size_t >::type parallel_type(parallel_typeSEXP);
+    Rcpp::traits::input_parameter< const IntegerVector& >::type parallel_arg(parallel_argSEXP);
+    rcpp_result_gen = Rcpp::wrap(gwr_generalized_fit(x, y, coords, family, bw, adaptive, kernel, longlat, p, theta, hatmatrix, intercept, optim_bw, optim_bw_criterion, parallel_type, parallel_arg));
+    return rcpp_result_gen;
+END_RCPP
+}
+// gwr_generalized_predict
+arma::mat gwr_generalized_predict(const arma::mat& pcoords, const arma::mat& x, const arma::vec& y, const arma::mat& coords, size_t family, double bw, bool adaptive, size_t kernel, bool longlat, double p, double theta, bool hatmatrix, bool intercept, size_t parallel_type, const IntegerVector& parallel_arg);
+RcppExport SEXP _GWmodel3_gwr_generalized_predict(SEXP pcoordsSEXP, SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP familySEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP hatmatrixSEXP, SEXP interceptSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type pcoords(pcoordsSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
+    Rcpp::traits::input_parameter< const arma::vec& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type coords(coordsSEXP);
+    Rcpp::traits::input_parameter< size_t >::type family(familySEXP);
+    Rcpp::traits::input_parameter< double >::type bw(bwSEXP);
+    Rcpp::traits::input_parameter< bool >::type adaptive(adaptiveSEXP);
+    Rcpp::traits::input_parameter< size_t >::type kernel(kernelSEXP);
+    Rcpp::traits::input_parameter< bool >::type longlat(longlatSEXP);
+    Rcpp::traits::input_parameter< double >::type p(pSEXP);
+    Rcpp::traits::input_parameter< double >::type theta(thetaSEXP);
+    Rcpp::traits::input_parameter< bool >::type hatmatrix(hatmatrixSEXP);
+    Rcpp::traits::input_parameter< bool >::type intercept(interceptSEXP);
+    Rcpp::traits::input_parameter< size_t >::type parallel_type(parallel_typeSEXP);
+    Rcpp::traits::input_parameter< const IntegerVector& >::type parallel_arg(parallel_argSEXP);
+    rcpp_result_gen = Rcpp::wrap(gwr_generalized_predict(pcoords, x, y, coords, family, bw, adaptive, kernel, longlat, p, theta, hatmatrix, intercept, parallel_type, parallel_arg));
+    return rcpp_result_gen;
+END_RCPP
+}
 // gwr_lcr_fit
 List gwr_lcr_fit(const arma::mat& x, const arma::vec& y, const arma::mat& coords, double bw, bool adaptive, size_t kernel, bool longlat, double p, double theta, double lambda, bool lambda_adjust, double cn_thresh, bool intercept, bool hatmatrix, size_t parallel_type, const IntegerVector& parallel_arg, bool optim_bw);
 RcppExport SEXP _GWmodel3_gwr_lcr_fit(SEXP xSEXP, SEXP ySEXP, SEXP coordsSEXP, SEXP bwSEXP, SEXP adaptiveSEXP, SEXP kernelSEXP, SEXP longlatSEXP, SEXP pSEXP, SEXP thetaSEXP, SEXP lambdaSEXP, SEXP lambda_adjustSEXP, SEXP cn_threshSEXP, SEXP interceptSEXP, SEXP hatmatrixSEXP, SEXP parallel_typeSEXP, SEXP parallel_argSEXP, SEXP optim_bwSEXP) {
@@ -231,6 +282,8 @@ static const R_CallMethodDef CallEntries[] = {
     {"_GWmodel3_gwpca_cal", (DL_FUNC) &_GWmodel3_gwpca_cal, 9},
     {"_GWmodel3_gwr_basic_fit", (DL_FUNC) &_GWmodel3_gwr_basic_fit, 22},
     {"_GWmodel3_gwr_basic_predict", (DL_FUNC) &_GWmodel3_gwr_basic_predict, 14},
+    {"_GWmodel3_gwr_generalized_fit", (DL_FUNC) &_GWmodel3_gwr_generalized_fit, 16},
+    {"_GWmodel3_gwr_generalized_predict", (DL_FUNC) &_GWmodel3_gwr_generalized_predict, 15},
     {"_GWmodel3_gwr_lcr_fit", (DL_FUNC) &_GWmodel3_gwr_lcr_fit, 17},
     {"_GWmodel3_gwr_multiscale_fit", (DL_FUNC) &_GWmodel3_gwr_multiscale_fit, 25},
     {NULL, NULL, 0}
diff --git a/src/gwr_generalized.cpp b/src/gwr_generalized.cpp
new file mode 100644
index 0000000..178da84
--- /dev/null
+++ b/src/gwr_generalized.cpp
@@ -0,0 +1,195 @@
+// [[Rcpp::depends(RcppArmadillo)]]
+#include <RcppArmadillo.h>
+#include "utils.h"
+#include "gwmodel.h"
+
+using namespace std;
+using namespace Rcpp;
+using namespace arma;
+using namespace gwm;
+
+// [[Rcpp::export]]
+List gwr_generalized_fit(
+    const arma::mat& x,
+    const arma::vec& y,
+    const arma::mat& coords,
+    size_t family,
+    double bw,
+    bool adaptive,
+    size_t kernel,
+    bool longlat,
+    double p,
+    double theta,
+    bool hatmatrix,
+    bool intercept,
+    bool optim_bw,
+    size_t optim_bw_criterion,
+    size_t parallel_type,
+    const IntegerVector& parallel_arg
+) {
+    vector<int> vpar_args = as< vector<int> >(IntegerVector(parallel_arg));
+
+    // Make Spatial Weight
+    BandwidthWeight bandwidth(bw, adaptive, BandwidthWeight::KernelFunctionType((size_t)kernel));
+    Distance* distance = nullptr;
+    if (longlat)
+    {
+        distance = new CRSDistance(true);
+    }
+    else
+    {
+        if (p == 2.0 && theta == 0.0)
+        {
+            distance = new CRSDistance(false);
+        }
+        else
+        {
+            distance = new MinkwoskiDistance(p, theta);
+        }
+    }
+    SpatialWeight spatial(&bandwidth, distance);
+
+    // Make Algorithm Object
+    GWRGeneralized algorithm;
+    algorithm.setCoords(coords);
+    algorithm.setDependentVariable(y);
+    algorithm.setIndependentVariables(x);
+    algorithm.setSpatialWeight(spatial);
+    algorithm.setHasHatMatrix(hatmatrix);
+    algorithm.setHasIntercept(intercept);
+    algorithm.setIsAutoselectBandwidth(optim_bw);
+    algorithm.setBandwidthSelectionCriterionType(GWRGeneralized::BandwidthSelectionCriterionType(optim_bw_criterion));
+
+    algorithm.setFamily(GWRGeneralized::Family(size_t(family)));
+
+    switch (ParallelType(size_t(parallel_type)))
+    {
+    case ParallelType::SerialOnly:
+        algorithm.setParallelType(ParallelType::SerialOnly);
+        break;
+#ifdef _OPENMP
+    case ParallelType::OpenMP:
+        algorithm.setParallelType(ParallelType::OpenMP);
+        algorithm.setOmpThreadNum(vpar_args[0]);
+        break;
+#endif
+    default:
+        algorithm.setParallelType(ParallelType::SerialOnly);
+        break;
+    }
+    try
+    {
+        algorithm.fit();
+    }
+    catch(const std::exception& e)
+    {
+        stop(e.what());
+    }
+
+    // Return Results
+    mat betas = algorithm.betas();
+    vec yhat = arma::sum(betas % x, 1);
+
+    List diagnostic =  List::create(
+        Named("RSS") = algorithm.getDiagnostic().RSS,
+        Named("AIC") = algorithm.getDiagnostic().AIC,
+        Named("AICc") = algorithm.getDiagnostic().AICc,
+        Named("RSquare") = algorithm.getDiagnostic().RSquare
+    );
+
+    List result_list = List::create(
+        Named("betas") = betas,
+        Named("yhat") = yhat,
+        Named("diagnostic") = diagnostic
+    );
+    if (optim_bw)
+    {
+        double bw_value = algorithm.spatialWeight().weight<BandwidthWeight>()->bandwidth();
+        result_list["bandwidth"] = wrap(bw_value);
+    }
+
+    return result_list;
+}
+
+
+// [[Rcpp::export]]
+arma::mat gwr_generalized_predict(
+    const arma::mat& pcoords,
+    const arma::mat& x,
+    const arma::vec& y,
+    const arma::mat& coords,
+    size_t family,
+    double bw,
+    bool adaptive,
+    size_t kernel,
+    bool longlat,
+    double p,
+    double theta,
+    bool hatmatrix,
+    bool intercept,
+    size_t parallel_type,
+    const IntegerVector& parallel_arg
+) {
+    // Convert data types
+    vector<int> vpar_args = as< vector<int> >(IntegerVector(parallel_arg));
+
+    // Make Spatial Weight
+    BandwidthWeight bandwidth(bw, adaptive, BandwidthWeight::KernelFunctionType((size_t)kernel));
+    Distance* distance = nullptr;
+    if (longlat)
+    {
+        distance = new CRSDistance(true);
+    }
+    else
+    {
+        if (p == 2.0 && theta == 0.0)
+        {
+            distance = new CRSDistance(false);
+        }
+        else
+        {
+            distance = new MinkwoskiDistance(p, theta);
+        }
+    }
+    SpatialWeight spatial(&bandwidth, distance);
+
+    // Make Algorithm Object
+    GWRGeneralized algorithm;
+    algorithm.setCoords(coords);
+    algorithm.setDependentVariable(y);
+    algorithm.setIndependentVariables(x);
+    algorithm.setSpatialWeight(spatial);
+    algorithm.setHasHatMatrix(hatmatrix);
+    algorithm.setHasIntercept(intercept);
+
+    algorithm.setFamily(GWRGeneralized::Family(size_t(family)));
+
+    switch (ParallelType(size_t(parallel_type)))
+    {
+    case ParallelType::SerialOnly:
+        algorithm.setParallelType(ParallelType::SerialOnly);
+        break;
+#ifdef _OPENMP
+    case ParallelType::OpenMP:
+        algorithm.setParallelType(ParallelType::OpenMP);
+        algorithm.setOmpThreadNum(vpar_args[0]);
+        break;
+#endif
+    default:
+        algorithm.setParallelType(ParallelType::SerialOnly);
+        break;
+    }
+
+    // Return Results
+    mat betas;
+    try
+    {
+        betas = algorithm.predict(pcoords);
+    }
+    catch(const std::exception& e)
+    {
+        stop(e.what());
+    }
+
+    return betas;
+}
\ No newline at end of file
diff --git a/tests/testthat/test-gwr_generalized.R b/tests/testthat/test-gwr_generalized.R
new file mode 100644
index 0000000..dcad187
--- /dev/null
+++ b/tests/testthat/test-gwr_generalized.R
@@ -0,0 +1,32 @@
+data(LondonHP)
+# m <- NULL
+
+test_that("GWR Generalized: poisson family", {
+  m1 <<- expect_no_error(
+    gwr_generalized(PURCHASE~FLOORSZ+UNEMPLOY, LondonHP, "poisson", 64, TRUE)
+  )
+})
+
+test_that("GWR Generalized: binomial family", {
+  expect_no_error(
+    gwr_generalized(TYPEDETCH~FLOORSZ+UNEMPLOY, LondonHP, "binomial", 64, TRUE)
+  )
+})
+
+test_that("GWR Generalized: Bandwidth selection", {
+  m2 <<- expect_no_error(
+    gwr_generalized(TYPEDETCH~FLOORSZ+UNEMPLOY, LondonHP, "binomial", bw="CV", TRUE)
+  )
+})
+
+test_that("GWR Generalized: Bandwidth selection, omp", {
+  expect_no_error(
+    gwr_generalized(TYPEDETCH~FLOORSZ+UNEMPLOY, LondonHP, "binomial", bw="AIC", TRUE, parallel_method = "omp", parallel_arg = 4)
+  )
+})
+
+# 测试问题可能在内核库的predict函数中，应该先设置一下mHasRegressionData，否则导致现在出现矩阵列数计算的错误。
+# 在内核库中测试通过predict后，再完善这部分内容。
+# test_that("GWR Generalized: predict", {
+#   predict(m1, LondonHP)
+# })
\ No newline at end of file