diff --git a/mllib/src/main/scala/org/apache/spark/ml/regression/GeneralizedLinearRegression.scala b/mllib/src/main/scala/org/apache/spark/ml/regression/GeneralizedLinearRegression.scala
new file mode 100644
index 0000000000000..68797d50ee3cc
--- /dev/null
+++ b/mllib/src/main/scala/org/apache/spark/ml/regression/GeneralizedLinearRegression.scala
@@ -0,0 +1,472 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *    http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.spark.ml.regression
+
+import breeze.stats.distributions.{Gaussian => GD}
+
+import org.apache.spark.Logging
+import org.apache.spark.annotation.{Experimental, Since}
+import org.apache.spark.ml.PredictorParams
+import org.apache.spark.ml.feature.Instance
+import org.apache.spark.ml.optim._
+import org.apache.spark.ml.param._
+import org.apache.spark.ml.param.shared._
+import org.apache.spark.ml.util.Identifiable
+import org.apache.spark.mllib.linalg.{BLAS, Vector}
+import org.apache.spark.rdd.RDD
+import org.apache.spark.sql.{DataFrame, Row}
+import org.apache.spark.sql.functions._
+
+/**
+ * Params for Generalized Linear Regression.
+ */
+private[regression] trait GeneralizedLinearRegressionParams extends PredictorParams
+  with HasFitIntercept with HasMaxIter with HasTol with HasRegParam with HasWeightCol
+  with HasSolver with Logging {
+
+  /**
+   * Param for the name of family which is a description of the error distribution
+   * to be used in the model.
+   * Supported options: "gaussian", "binomial", "poisson" and "gamma".
+   * @group param
+   */
+  @Since("2.0.0")
+  final val family: Param[String] = new Param(this, "family",
+    "the name of family which is a description of the error distribution to be used in the model",
+    ParamValidators.inArray[String](GeneralizedLinearRegression.supportedFamilies.toArray))
+
+  /** @group getParam */
+  @Since("2.0.0")
+  def getFamily: String = $(family)
+
+  /**
+   * Param for the name of the model link function.
+   * Supported options: "identity", "log", "inverse", "logit", "probit", "cloglog" and "sqrt".
+   * @group param
+   */
+  @Since("2.0.0")
+  final val link: Param[String] = new Param(this, "link", "the name of the model link function",
+    ParamValidators.inArray[String](GeneralizedLinearRegression.supportedLinks.toArray))
+
+  /** @group getParam */
+  @Since("2.0.0")
+  def getLink: String = $(link)
+
+  @Since("2.0.0")
+  override def validateParams(): Unit = {
+    require(GeneralizedLinearRegression.supportedFamilyLinkPairs.contains($(family) -> $(link)),
+      s"Generalized Linear Regression with ${$(family)} family does not support ${$(link)} " +
+        s"link function.")
+  }
+}
+
+/**
+ * :: Experimental ::
+ *
+ * Fit a Generalized Linear Model ([[https://en.wikipedia.org/wiki/Generalized_linear_model]])
+ * specified by giving a symbolic description of the linear predictor and
+ * a description of the error distribution.
+ */
+@Experimental
+@Since("2.0.0")
+class GeneralizedLinearRegression @Since("2.0.0") (@Since("2.0.0") override val uid: String)
+  extends Regressor[Vector, GeneralizedLinearRegression, GeneralizedLinearRegressionModel]
+  with GeneralizedLinearRegressionParams with Logging {
+
+  @Since("2.0.0")
+  def this() = this(Identifiable.randomUID("genLinReg"))
+
+  /**
+   * Set the name of family which is a description of the error distribution
+   * to be used in the model.
+   * @group setParam
+   */
+  @Since("2.0.0")
+  def setFamily(value: String): this.type = set(family, value)
+
+  /**
+   * Set the name of the model link function.
+   * @group setParam
+   */
+  @Since("2.0.0")
+  def setLink(value: String): this.type = set(link, value)
+
+  /**
+   * Set if we should fit the intercept.
+   * Default is true.
+   * @group setParam
+   */
+  @Since("2.0.0")
+  def setFitIntercept(value: Boolean): this.type = set(fitIntercept, value)
+  setDefault(fitIntercept -> true)
+
+  /**
+   * Set the maximum number of iterations.
+   * Default is 100.
+   * @group setParam
+   */
+  @Since("2.0.0")
+  def setMaxIter(value: Int): this.type = set(maxIter, value)
+  setDefault(maxIter -> 100)
+
+  /**
+   * Set the convergence tolerance of iterations.
+   * Smaller value will lead to higher accuracy with the cost of more iterations.
+   * Default is 1E-6.
+   * @group setParam
+   */
+  @Since("2.0.0")
+  def setTol(value: Double): this.type = set(tol, value)
+  setDefault(tol -> 1E-6)
+
+  /**
+   * Set the regularization parameter.
+   * Default is 0.0.
+   * @group setParam
+   */
+  @Since("2.0.0")
+  def setRegParam(value: Double): this.type = set(regParam, value)
+  setDefault(regParam -> 0.0)
+
+  /**
+   * Whether to over-/under-sample training instances according to the given weights in weightCol.
+   * If empty, all instances are treated equally (weight 1.0).
+   * Default is empty, so all instances have weight one.
+   * @group setParam
+   */
+  @Since("2.0.0")
+  def setWeightCol(value: String): this.type = set(weightCol, value)
+  setDefault(weightCol -> "")
+
+  /**
+   * Set the solver algorithm used for optimization.
+   * Currently only support "irls" which is also the default solver.
+   * @group setParam
+   */
+  @Since("2.0.0")
+  def setSolver(value: String): this.type = set(solver, value)
+  setDefault(solver -> "irls")
+
+  override protected def train(dataset: DataFrame): GeneralizedLinearRegressionModel = {
+    val familyLink = $(family) match {
+      case "gaussian" => Gaussian($(link))
+      case "binomial" => Binomial($(link))
+      case "poisson" => Poisson($(link))
+      case "gamma" => Gamma($(link))
+    }
+
+    val w = if ($(weightCol).isEmpty) lit(1.0) else col($(weightCol))
+    val instances: RDD[Instance] = dataset.select(
+      col($(labelCol)), w, col($(featuresCol))).map {
+      case Row(label: Double, weight: Double, features: Vector) =>
+        Instance(label, weight, features)
+    }
+
+    if ($(family) == "gaussian" && $(link) == "identity") {
+      val optimizer = new WeightedLeastSquares($(fitIntercept), $(regParam),
+        standardizeFeatures = true, standardizeLabel = true)
+      val wlsModel = optimizer.fit(instances)
+      val model = copyValues(new GeneralizedLinearRegressionModel(uid,
+        wlsModel.coefficients, wlsModel.intercept).setParent(this))
+      return model
+    }
+
+    val newInstances = instances.map { instance =>
+      val mu = familyLink.initialize(instance.label, instance.weight)
+      val eta = familyLink.predict(mu)
+      Instance(eta, instance.weight, instance.features)
+    }
+
+    val initialModel = new WeightedLeastSquares($(fitIntercept), $(regParam),
+      standardizeFeatures = true, standardizeLabel = true).fit(newInstances)
+
+    val reweightFunc: (Instance, WeightedLeastSquaresModel) => (Double, Double) = {
+      (instance: Instance, model: WeightedLeastSquaresModel) => {
+        val eta = model.predict(instance.features)
+        val mu = familyLink.fitted(eta)
+        val z = familyLink.adjusted(instance.label, mu, eta)
+        val w = familyLink.weights(mu) * instance.weight
+        (z, w)
+      }
+    }
+
+    val optimizer = new IterativelyReweightedLeastSquares(initialModel, reweightFunc,
+      $(fitIntercept), $(regParam), $(maxIter), $(tol))
+
+    val irlsModel = optimizer.fit(instances)
+
+    val model = copyValues(new GeneralizedLinearRegressionModel(uid,
+      irlsModel.coefficients, irlsModel.intercept))
+    model
+  }
+
+  @Since("2.0.0")
+  override def copy(extra: ParamMap): GeneralizedLinearRegression = defaultCopy(extra)
+}
+
+@Since("2.0.0")
+object GeneralizedLinearRegression {
+
+  /** Set of families that GeneralizedLinearRegression supports */
+  private[ml] val supportedFamilies = Set("gaussian", "binomial", "poisson", "gamma")
+
+  /** Set of links that GeneralizedLinearRegression supports */
+  private[ml] val supportedLinks = Set("identity", "log", "inverse", "logit", "probit",
+    "cloglog", "sqrt")
+
+  /** Set of family and link pairs that GeneralizedLinearRegression supports */
+  private[ml] val supportedFamilyLinkPairs = Set(
+    "gaussian" -> "identity", "gaussian" -> "log", "gaussian" -> "inverse",
+    "binomial" -> "logit", "binomial" -> "probit", "binomial" -> "cloglog",
+    "poisson" -> "log", "poisson" -> "identity", "poisson" -> "sqrt",
+    "gamma" -> "inverse", "gamma" -> "identity", "gamma" -> "log"
+  )
+}
+
+/**
+ * :: Experimental ::
+ * Model produced by [[GeneralizedLinearRegression]].
+ */
+@Experimental
+@Since("2.0.0")
+class GeneralizedLinearRegressionModel private[ml] (
+    @Since("2.0.0") override val uid: String,
+    @Since("2.0.0") val coefficients: Vector,
+    @Since("2.0.0") val intercept: Double)
+  extends RegressionModel[Vector, GeneralizedLinearRegressionModel]
+  with GeneralizedLinearRegressionParams {
+
+  private lazy val familyLink = $(family) match {
+    case "gaussian" => Gaussian($(link))
+    case "binomial" => Binomial($(link))
+    case "poisson" => Poisson($(link))
+    case "gamma" => Gamma($(link))
+  }
+
+  override protected def predict(features: Vector): Double = {
+    val eta = BLAS.dot(features, coefficients) + intercept
+    familyLink.fitted(eta)
+  }
+
+  @Since("2.0.0")
+  override def copy(extra: ParamMap): GeneralizedLinearRegressionModel = {
+    copyValues(new GeneralizedLinearRegressionModel(uid, coefficients, intercept), extra)
+      .setParent(parent)
+  }
+}
+
+/**
+ * A description of the error distribution and link function to be used in the model.
+ * @param link a link function instance
+ */
+private[ml] abstract class Family(val link: Link) extends Serializable {
+
+  /** Initialize the starting value for mu. */
+  def initialize(y: Double, weight: Double): Double
+
+  /** The variance of mu to its mean. */
+  def variance(mu: Double): Double
+
+  /** Weights for IRLS steps. */
+  def weights(mu: Double): Double = {
+    1.0 / (math.pow(this.link.deriv(mu), 2.0) * this.variance(mu))
+  }
+
+  /** The adjusted response variable. */
+  def adjusted(y: Double, mu: Double, eta: Double): Double = {
+    eta + (y - mu) * link.deriv(mu)
+  }
+
+  /** Linear predictors based on given mu. */
+  def predict(mu: Double): Double = this.link.link(mu)
+
+  /** Fitted values based on linear predictors eta. */
+  def fitted(eta: Double): Double = this.link.unlink(eta)
+}
+
+/**
+ * Gaussian exponential family distribution.
+ * The default link for the Gaussian family is the identity link.
+ * @param link a link function instance
+ */
+private[ml] class Gaussian(link: Link = Identity) extends Family(link) {
+
+  override def initialize(y: Double, weight: Double): Double = y
+
+  def variance(mu: Double): Double = 1.0
+}
+
+private[ml] object Gaussian {
+
+  def apply(link: String): Gaussian = {
+    link match {
+      case "identity" => new Gaussian(Identity)
+      case "log" => new Gaussian(Log)
+      case "inverse" => new Gaussian(Inverse)
+    }
+  }
+}
+
+/**
+ * Binomial exponential family distribution.
+ * The default link for the Binomial family is the logit link.
+ * @param link a link function instance
+ */
+private[ml] class Binomial(link: Link = Logit) extends Family(link) {
+
+  override def initialize(y: Double, weight: Double): Double = {
+    (weight * y + 0.5) / (weight + 1.0)
+  }
+
+  override def variance(mu: Double): Double = mu * (1 - mu)
+}
+
+private[ml] object Binomial {
+
+  def apply(link: String): Binomial = {
+    link match {
+      case "logit" => new Binomial(Logit)
+      case "probit" => new Binomial(Probit)
+      case "cloglog" => new Binomial(CLogLog)
+    }
+  }
+}
+
+/**
+ * Poisson exponential family distribution.
+ * The default link for the Poisson family is the log link.
+ * @param link a link function instance
+ */
+private[ml] class Poisson(link: Link = Log) extends Family(link) {
+
+  override def initialize(y: Double, weight: Double): Double = y + 0.1
+
+  override def variance(mu: Double): Double = mu
+}
+
+private[ml] object Poisson {
+
+  def apply(link: String): Poisson = {
+    link match {
+      case "log" => new Poisson(Log)
+      case "sqrt" => new Poisson(Sqrt)
+      case "identity" => new Poisson(Identity)
+    }
+  }
+}
+
+/**
+ * Gamma exponential family distribution.
+ * The default link for the Gamma family is the inverse link.
+ * @param link a link function instance
+ */
+private[ml] class Gamma(link: Link = Log) extends Family(link) {
+
+  override def initialize(y: Double, weight: Double): Double = y
+
+  override def variance(mu: Double): Double = math.pow(mu, 2.0)
+}
+
+private[ml] object Gamma {
+
+  def apply(link: String): Gamma = {
+    link match {
+      case "inverse" => new Gamma(Inverse)
+      case "identity" => new Gamma(Identity)
+      case "log" => new Gamma(Log)
+    }
+  }
+}
+
+/**
+ * A description of the link function to be used in the model.
+ */
+private[ml] trait Link extends Serializable {
+
+  /** The link function. */
+  def link(mu: Double): Double
+
+  /** The derivative function. */
+  def deriv(mu: Double): Double
+
+  /** The inverse link function. */
+  def unlink(eta: Double): Double
+}
+
+private[ml] object Identity extends Link {
+
+  override def link(mu: Double): Double = mu
+
+  override def deriv(mu: Double): Double = 1.0
+
+  override def unlink(eta: Double): Double = eta
+}
+
+private[ml] object Logit extends Link {
+
+  override def link(mu: Double): Double = math.log(mu / (1.0 - mu))
+
+  override def deriv(mu: Double): Double = 1.0 / (mu * (1.0 - mu))
+
+  override def unlink(eta: Double): Double = 1.0 / (1.0 + math.exp(-1.0 * eta))
+}
+
+private[ml] object Log extends Link {
+
+  override def link(mu: Double): Double = math.log(mu)
+
+  override def deriv(mu: Double): Double = 1.0 / mu
+
+  override def unlink(eta: Double): Double = math.exp(eta)
+}
+
+private[ml] object Inverse extends Link {
+
+  override def link(mu: Double): Double = 1.0 / mu
+
+  override def deriv(mu: Double): Double = -1.0 * math.pow(mu, -2.0)
+
+  override def unlink(eta: Double): Double = 1.0 / eta
+}
+
+private[ml] object Probit extends Link {
+
+  override def link(mu: Double): Double = GD(0.0, 1.0).icdf(mu)
+
+  override def deriv(mu: Double): Double = 1.0 / GD(0.0, 1.0).pdf(GD(0.0, 1.0).icdf(mu))
+
+  override def unlink(eta: Double): Double = GD(0.0, 1.0).cdf(eta)
+}
+
+private[ml] object CLogLog extends Link {
+
+  override def link(mu: Double): Double = math.log(-1.0 * math.log(1 - mu))
+
+  override def deriv(mu: Double): Double = 1.0 / ((mu - 1.0) * math.log(1.0 - mu))
+
+  override def unlink(eta: Double): Double = 1.0 - math.exp(-1.0 * math.exp(eta))
+}
+
+private[ml] object Sqrt extends Link {
+
+  override def link(mu: Double): Double = math.sqrt(mu)
+
+  override def deriv(mu: Double): Double = 1.0 / (2.0 * math.sqrt(mu))
+
+  override def unlink(eta: Double): Double = math.pow(eta, 2.0)
+}
diff --git a/mllib/src/test/scala/org/apache/spark/ml/regression/GeneralizedLinearRegressionSuite.scala b/mllib/src/test/scala/org/apache/spark/ml/regression/GeneralizedLinearRegressionSuite.scala
new file mode 100644
index 0000000000000..abfaad57d7ba8
--- /dev/null
+++ b/mllib/src/test/scala/org/apache/spark/ml/regression/GeneralizedLinearRegressionSuite.scala
@@ -0,0 +1,448 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *    http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.spark.ml.regression
+
+import scala.util.Random
+
+import org.apache.spark.SparkFunSuite
+import org.apache.spark.ml.param.ParamsSuite
+import org.apache.spark.ml.util.MLTestingUtils
+import org.apache.spark.mllib.classification.LogisticRegressionSuite._
+import org.apache.spark.mllib.linalg.{BLAS, DenseVector, Vectors}
+import org.apache.spark.mllib.random._
+import org.apache.spark.mllib.regression.LabeledPoint
+import org.apache.spark.mllib.util.MLlibTestSparkContext
+import org.apache.spark.mllib.util.TestingUtils._
+import org.apache.spark.sql.{DataFrame, Row}
+
+class GeneralizedLinearRegressionSuite extends SparkFunSuite with MLlibTestSparkContext {
+
+  private val seed: Int = 42
+  @transient var datasetGaussianIdentity: DataFrame = _
+  @transient var datasetGaussianLog: DataFrame = _
+  @transient var datasetGaussianInverse: DataFrame = _
+  @transient var datasetBinomial: DataFrame = _
+  @transient var datasetPoissonLog: DataFrame = _
+  @transient var datasetPoissonIdentity: DataFrame = _
+  @transient var datasetPoissonSqrt: DataFrame = _
+  @transient var datasetGammaInverse: DataFrame = _
+  @transient var datasetGammaIdentity: DataFrame = _
+  @transient var datasetGammaLog: DataFrame = _
+
+  override def beforeAll(): Unit = {
+    super.beforeAll()
+
+    import GeneralizedLinearRegressionSuite._
+
+    datasetGaussianIdentity = sqlContext.createDataFrame(
+      sc.parallelize(generateGeneralizedLinearRegressionInput(
+        intercept = 2.5, coefficients = Array(2.2, 0.6), xMean = Array(2.9, 10.5),
+        xVariance = Array(0.7, 1.2), nPoints = 10000, seed, eps = 0.01,
+        family = "gaussian", link = "identity"), 2))
+
+    datasetGaussianLog = sqlContext.createDataFrame(
+      sc.parallelize(generateGeneralizedLinearRegressionInput(
+        intercept = 0.25, coefficients = Array(0.22, 0.06), xMean = Array(2.9, 10.5),
+        xVariance = Array(0.7, 1.2), nPoints = 10000, seed, eps = 0.01,
+        family = "gaussian", link = "log"), 2))
+
+    datasetGaussianInverse = sqlContext.createDataFrame(
+      sc.parallelize(generateGeneralizedLinearRegressionInput(
+        intercept = 2.5, coefficients = Array(2.2, 0.6), xMean = Array(2.9, 10.5),
+        xVariance = Array(0.7, 1.2), nPoints = 10000, seed, eps = 0.01,
+        family = "gaussian", link = "inverse"), 2))
+
+    datasetBinomial = {
+      val nPoints = 10000
+      val coefficients = Array(-0.57997, 0.912083, -0.371077, -0.819866, 2.688191)
+      val xMean = Array(5.843, 3.057, 3.758, 1.199)
+      val xVariance = Array(0.6856, 0.1899, 3.116, 0.581)
+
+      val testData =
+        generateMultinomialLogisticInput(coefficients, xMean, xVariance, true, nPoints, seed)
+
+      sqlContext.createDataFrame(sc.parallelize(testData, 4))
+    }
+
+    datasetPoissonLog = sqlContext.createDataFrame(
+      sc.parallelize(generateGeneralizedLinearRegressionInput(
+        intercept = 0.25, coefficients = Array(0.22, 0.06), xMean = Array(2.9, 10.5),
+        xVariance = Array(0.7, 1.2), nPoints = 10000, seed, eps = 0.01,
+        family = "poisson", link = "log"), 2))
+
+    datasetPoissonIdentity = sqlContext.createDataFrame(
+      sc.parallelize(generateGeneralizedLinearRegressionInput(
+        intercept = 2.5, coefficients = Array(2.2, 0.6), xMean = Array(2.9, 10.5),
+        xVariance = Array(0.7, 1.2), nPoints = 10000, seed, eps = 0.01,
+        family = "poisson", link = "identity"), 2))
+
+    datasetPoissonSqrt = sqlContext.createDataFrame(
+      sc.parallelize(generateGeneralizedLinearRegressionInput(
+        intercept = 2.5, coefficients = Array(2.2, 0.6), xMean = Array(2.9, 10.5),
+        xVariance = Array(0.7, 1.2), nPoints = 10000, seed, eps = 0.01,
+        family = "poisson", link = "sqrt"), 2))
+
+    datasetGammaInverse = sqlContext.createDataFrame(
+      sc.parallelize(generateGeneralizedLinearRegressionInput(
+        intercept = 2.5, coefficients = Array(2.2, 0.6), xMean = Array(2.9, 10.5),
+        xVariance = Array(0.7, 1.2), nPoints = 10000, seed, eps = 0.01,
+        family = "gamma", link = "inverse"), 2))
+
+    datasetGammaIdentity = sqlContext.createDataFrame(
+      sc.parallelize(generateGeneralizedLinearRegressionInput(
+        intercept = 2.5, coefficients = Array(2.2, 0.6), xMean = Array(2.9, 10.5),
+        xVariance = Array(0.7, 1.2), nPoints = 10000, seed, eps = 0.01,
+        family = "gamma", link = "identity"), 2))
+
+    datasetGammaLog = sqlContext.createDataFrame(
+      sc.parallelize(generateGeneralizedLinearRegressionInput(
+        intercept = 0.25, coefficients = Array(0.22, 0.06), xMean = Array(2.9, 10.5),
+        xVariance = Array(0.7, 1.2), nPoints = 10000, seed, eps = 0.01,
+        family = "gamma", link = "log"), 2))
+  }
+
+  test("params") {
+    ParamsSuite.checkParams(new GeneralizedLinearRegression)
+    val model = new GeneralizedLinearRegressionModel("genLinReg", Vectors.dense(0.0), 0.0)
+    ParamsSuite.checkParams(model)
+  }
+
+  test("generalized linear regression: default params") {
+    val glr = new GeneralizedLinearRegression
+    assert(glr.getLabelCol === "label")
+    assert(glr.getFeaturesCol === "features")
+    assert(glr.getPredictionCol === "prediction")
+    assert(glr.getFitIntercept)
+    assert(glr.getMaxIter === 100)
+    assert(glr.getTol === 1E-6)
+    assert(glr.getWeightCol === "")
+    assert(glr.getRegParam === 0.0)
+    assert(glr.getSolver == "irls")
+    val model = glr.setFamily("gaussian").setLink("identity")
+      .fit(datasetGaussianIdentity)
+
+    // copied model must have the same parent.
+    MLTestingUtils.checkCopy(model)
+    model.transform(datasetGaussianIdentity)
+      .select("label", "prediction")
+      .collect()
+
+    assert(model.getFeaturesCol === "features")
+    assert(model.getPredictionCol === "prediction")
+    assert(model.intercept !== 0.0)
+    assert(model.hasParent)
+    assert(model.getFamily === "gaussian")
+    assert(model.getLink === "identity")
+  }
+
+  test("generalized linear regression: gaussian family") {
+    /*
+       R code:
+       f1 <- data$V1 ~ data$V2 + data$V3 - 1
+       f2 <- data$V1 ~ data$V2 + data$V3
+
+       data <- read.csv("path", header=FALSE)
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family="gaussian", data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] 2.2960999 0.8087933
+       [1] 2.5002642 2.2000403 0.5999485
+
+       data <- read.csv("path", header=FALSE)
+       model1 <- glm(f1, family=gaussian(link=log), data=data, start=c(0,0))
+       model2 <- glm(f2, family=gaussian(link=log), data=data, start=c(0,0,0))
+       print(as.vector(coef(model1)))
+       print(as.vector(coef(model2)))
+
+       [1] 0.23069326 0.07993778
+       [1] 0.25001858 0.22002452 0.05998789
+
+       data <- read.csv("path", header=FALSE)
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family=gaussian(link=inverse), data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] 2.3010179 0.8198976
+       [1] 2.4108902 2.2130248 0.6086152
+     */
+
+    val expected = Seq(
+      Vectors.dense(0.0, 2.2960999, 0.8087933),
+      Vectors.dense(2.5002642, 2.2000403, 0.5999485),
+      Vectors.dense(0.0, 0.23069326, 0.07993778),
+      Vectors.dense(0.25001858, 0.22002452, 0.05998789),
+      Vectors.dense(0.0, 2.3010179, 0.8198976),
+      Vectors.dense(2.4108902, 2.2130248, 0.6086152))
+
+    var idx = 0
+    for ((link, dataset) <- Seq(("identity", datasetGaussianIdentity), ("log", datasetGaussianLog),
+      ("inverse", datasetGaussianInverse))) {
+      for (fitIntercept <- Seq(false, true)) {
+        val trainer = new GeneralizedLinearRegression().setFamily("gaussian").setLink(link)
+          .setFitIntercept(fitIntercept)
+        val model = trainer.fit(dataset)
+        val actual = Vectors.dense(model.intercept, model.coefficients(0), model.coefficients(1))
+        assert(actual ~== expected(idx) absTol 1e-4)
+
+        val familyLink = Gaussian(link)
+        model.transform(dataset).select("features", "prediction").collect().foreach {
+          case Row(features: DenseVector, prediction1: Double) =>
+            val eta = BLAS.dot(features, model.coefficients) + model.intercept
+            val prediction2 = familyLink.fitted(eta)
+            assert(prediction1 ~== prediction2 relTol 1E-5)
+        }
+
+        idx += 1
+      }
+    }
+  }
+
+  test("generalized linear regression: binomial family") {
+    /*
+       R code:
+       f1 <- data$V1 ~ data$V2 + data$V3 + data$V4 + data$V5 - 1
+       f2 <- data$V1 ~ data$V2 + data$V3 + data$V4 + data$V5
+       data <- read.csv("path", header=FALSE)
+
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family="binomial", data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] -0.3560284  1.3010002 -0.3570805 -0.7406762
+       [1]  2.8367406 -0.5896187  0.8931655 -0.3925169 -0.7996989
+
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family=binomial(link=probit), data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] -0.2134390  0.7800646 -0.2144267 -0.4438358
+       [1]  1.6995366 -0.3524694  0.5332651 -0.2352985 -0.4780850
+
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family=binomial(link=cloglog), data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] -0.2832198  0.8434144 -0.2524727 -0.5293452
+       [1]  1.5063590 -0.4038015  0.6133664 -0.2687882 -0.5541758
+     */
+    val expected = Seq(
+      Vectors.dense(0.0, -0.3560284, 1.3010002, -0.3570805, -0.7406762),
+      Vectors.dense(2.8367406, -0.5896187, 0.8931655, -0.3925169, -0.7996989),
+      Vectors.dense(0.0, -0.2134390, 0.7800646, -0.2144267, -0.4438358),
+      Vectors.dense(1.6995366, -0.3524694, 0.5332651, -0.2352985, -0.4780850),
+      Vectors.dense(0.0, -0.2832198, 0.8434144, -0.2524727, -0.5293452),
+      Vectors.dense(1.5063590, -0.4038015, 0.6133664, -0.2687882, -0.5541758))
+
+    var idx = 0
+    for ((link, dataset) <- Seq(("logit", datasetBinomial), ("probit", datasetBinomial),
+      ("cloglog", datasetBinomial))) {
+      for (fitIntercept <- Seq(false, true)) {
+        val trainer = new GeneralizedLinearRegression().setFamily("binomial").setLink(link)
+          .setFitIntercept(fitIntercept)
+        val model = trainer.fit(dataset)
+        val actual = Vectors.dense(model.intercept, model.coefficients(0), model.coefficients(1),
+          model.coefficients(2), model.coefficients(3))
+        assert(actual ~== expected(idx) absTol 1e-4)
+
+        val familyLink = Binomial(link)
+        model.transform(dataset).select("features", "prediction").collect().foreach {
+          case Row(features: DenseVector, prediction1: Double) =>
+            val eta = BLAS.dot(features, model.coefficients) + model.intercept
+            val prediction2 = familyLink.fitted(eta)
+            assert(prediction1 ~== prediction2 relTol 1E-5)
+        }
+
+        idx += 1
+      }
+    }
+  }
+
+  test("generalized linear regression: poisson family") {
+    /*
+       R code:
+       f1 <- data$V1 ~ data$V2 + data$V3 - 1
+       f2 <- data$V1 ~ data$V2 + data$V3
+
+       data <- read.csv("path", header=FALSE)
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family="poisson", data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] 0.22999393 0.08047088
+       [1] 0.25022353 0.21998599 0.05998621
+
+       data <- read.csv("path", header=FALSE)
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family=poisson(link=identity), data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] 2.2929501 0.8119415
+       [1] 2.5012730 2.1999407 0.5999107
+
+       data <- read.csv("path", header=FALSE)
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family=poisson(link=sqrt), data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] 2.2958947 0.8090515
+       [1] 2.5000480 2.1999972 0.5999968
+     */
+    val expected = Seq(
+      Vectors.dense(0.0, 0.22999393, 0.08047088),
+      Vectors.dense(0.25022353, 0.21998599, 0.05998621),
+      Vectors.dense(0.0, 2.2929501, 0.8119415),
+      Vectors.dense(2.5012730, 2.1999407, 0.5999107),
+      Vectors.dense(0.0, 2.2958947, 0.8090515),
+      Vectors.dense(2.5000480, 2.1999972, 0.5999968))
+
+    var idx = 0
+    for ((link, dataset) <- Seq(("log", datasetPoissonLog), ("identity", datasetPoissonIdentity),
+      ("sqrt", datasetPoissonSqrt))) {
+      for (fitIntercept <- Seq(false, true)) {
+        val trainer = new GeneralizedLinearRegression().setFamily("poisson").setLink(link)
+          .setFitIntercept(fitIntercept)
+        val model = trainer.fit(dataset)
+        val actual = Vectors.dense(model.intercept, model.coefficients(0), model.coefficients(1))
+        assert(actual ~== expected(idx) absTol 1e-4)
+
+        val familyLink = Poisson(link)
+        model.transform(dataset).select("features", "prediction").collect().foreach {
+          case Row(features: DenseVector, prediction1: Double) =>
+            val eta = BLAS.dot(features, model.coefficients) + model.intercept
+            val prediction2 = familyLink.fitted(eta)
+            assert(prediction1 ~== prediction2 relTol 1E-5)
+        }
+
+        idx += 1
+      }
+    }
+  }
+
+  test("generalized linear regression: gamma family") {
+    /*
+       R code:
+       f1 <- data$V1 ~ data$V2 + data$V3 - 1
+       f2 <- data$V1 ~ data$V2 + data$V3
+
+       data <- read.csv("path", header=FALSE)
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family="Gamma", data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] 2.3392419 0.8058058
+       [1] 2.3507700 2.2533574 0.6042991
+
+       data <- read.csv("path", header=FALSE)
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family=Gamma(link=identity), data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] 2.2908883 0.8147796
+       [1] 2.5002406 2.1998346 0.6000059
+
+       data <- read.csv("path", header=FALSE)
+       for (formula in c(f1, f2)) {
+         model <- glm(formula, family=Gamma(link=log), data=data)
+         print(as.vector(coef(model)))
+       }
+
+       [1] 0.22958970 0.08091066
+       [1] 0.25003210 0.21996957 0.06000215
+     */
+    val expected = Seq(
+      Vectors.dense(0.0, 2.3392419, 0.8058058),
+      Vectors.dense(2.3507700, 2.2533574, 0.6042991),
+      Vectors.dense(0.0, 2.2908883, 0.8147796),
+      Vectors.dense(2.5002406, 2.1998346, 0.6000059),
+      Vectors.dense(0.0, 0.22958970, 0.08091066),
+      Vectors.dense(0.25003210, 0.21996957, 0.06000215))
+
+    var idx = 0
+    for ((link, dataset) <- Seq(("inverse", datasetGammaInverse),
+      ("identity", datasetGammaIdentity), ("log", datasetGammaLog))) {
+      for (fitIntercept <- Seq(false, true)) {
+        val trainer = new GeneralizedLinearRegression().setFamily("gamma").setLink(link)
+          .setFitIntercept(fitIntercept)
+        val model = trainer.fit(dataset)
+        val actual = Vectors.dense(model.intercept, model.coefficients(0), model.coefficients(1))
+        assert(actual ~== expected(idx) absTol 1e-4)
+
+        val familyLink = Gamma(link)
+        model.transform(dataset).select("features", "prediction").collect().foreach {
+          case Row(features: DenseVector, prediction1: Double) =>
+            val eta = BLAS.dot(features, model.coefficients) + model.intercept
+            val prediction2 = familyLink.fitted(eta)
+            assert(prediction1 ~== prediction2 relTol 1E-5)
+        }
+
+        idx += 1
+      }
+    }
+  }
+}
+
+object GeneralizedLinearRegressionSuite {
+
+  def generateGeneralizedLinearRegressionInput(
+      intercept: Double,
+      coefficients: Array[Double],
+      xMean: Array[Double],
+      xVariance: Array[Double],
+      nPoints: Int,
+      seed: Int,
+      eps: Double,
+      family: String,
+      link: String): Seq[LabeledPoint] = {
+
+    val rnd = new Random(seed)
+    def rndElement(i: Int) = {
+      (rnd.nextDouble() - 0.5) * math.sqrt(12.0 * xVariance(i)) + xMean(i)
+    }
+    val (generator, mean) = family match {
+      case "gaussian" => (new StandardNormalGenerator, 0.0)
+      case "poisson" => (new PoissonGenerator(1.0), 1.0)
+      case "gamma" => (new GammaGenerator(1.0, 1.0), 1.0)
+    }
+    generator.setSeed(seed)
+
+    (0 until nPoints).map { _ =>
+      val features = Vectors.dense(coefficients.indices.map { rndElement(_) }.toArray)
+      val eta = BLAS.dot(Vectors.dense(coefficients), features) + intercept
+      val mu = link match {
+        case "identity" => eta
+        case "log" => math.exp(eta)
+        case "sqrt" => math.pow(eta, 2.0)
+        case "inverse" => 1.0 / eta
+      }
+      val label = mu + eps * (generator.nextValue() - mean)
+      // Return LabeledPoints with DenseVector
+      LabeledPoint(label, features)
+    }
+  }
+}